Fully quantum description of the Zundel ion: combining variational quantum Monte Carlo with path integral Langevin dynamics

TL;DR

This paper presents a new fully quantum computational approach combining variational quantum Monte Carlo and path integral Langevin dynamics to accurately simulate coupled electron-ion systems, validated on the Zundel ion.

Contribution

It introduces a novel algorithm that integrates variational QMC with path integral Langevin dynamics, enabling unbiased, efficient quantum nuclear simulations with controlled statistical errors.

Findings

Nuclear quantum effects dominate over thermal effects in the Zundel ion.

Quantum tunneling increases proton mobility beyond classical predictions.

The method accurately reproduces results comparable to high-level quantum chemistry approaches.

Abstract

We introduce a novel approach for a fully quantum description of coupled electron-ion systems from first principles. It combines the variational quantum Monte Carlo (QMC) solution of the electronic part with the path integral (PI) formalism for the quantum nuclear dynamics. On the one hand, the PI molecular dynamics includes nuclear quantum effects by adding a set of fictitious classical particles (beads) aimed at reproducing nuclear quantum fluctuations via a harmonic kinetic term. On the other hand, variational QMC can provide Born-Oppenheimer (BO) potential energy surfaces with a precision comparable to the most advanced post Hartree-Fock approaches, and with a favorable scaling with the system size. To deal with the intrinsic QMC noise, we generalize the PI molecular dynamics using a Langevin thermostat correlated according to the covariance matrix of QMC nuclear forces. The variational parameters of the QMC wave function are evolved during the nuclear dynamics, such that the BO potential energy surface is unbiased. Statistical errors on the wave function parameters are reduced by resorting to bead grouping average, which we show to be accurate and well controlled. Our general algorithm relies on a Trotter breakup between the dynamics driven by ionic forces and the one set by the harmonic interbead couplings. The latter is exactly integrated even in presence of the Langevin thermostat, thanks to the mapping onto an Ornstein-Uhlenbeck process. This framework turns out to be very efficient also in the case of deterministic ionic forces. The new implementation is validated on the Zundel ion by direct comparison with standard PI Langevin dynamics calculations made with a coupled cluster potential energy surface. Nuclear quantum effects are confirmed to be dominant over thermal effects well beyond room temperature giving the excess proton an increased mobility by quantum tunneling.

Equations

$$H=\sum _{i=1}^{3N}\frac{p_i^2}{2m_i}+V\left(q_1,...,q_{3N}\right),$$

$$q_i$$

$$p_i$$

$${\eta }_i$$

$$f_i$$

$$\dot{\mathbf{p}}$$

$$\dot{\mathbf{q}}$$

$$\mathbf{\gamma }=\frac{\mathbf{\alpha }}{2k_{B}T},$$

$${\alpha }_{ij}\delta (t-t^{\prime })=⟨{\eta }_i(t){\eta }_j(t^{\prime })⟩$$

$$\mathbf{X}=\left(\begin{array}{c}\mathbf{p}\\ \mathbf{q}\end{array}\right),$$

$$\dot{\mathbf{X}}=-\hat{\mathbf{\gamma }}\mathbf{X}+\mathbf{F}+\mathbf{E},$$

$$\hat{\mathbf{\gamma }}=\left(\begin{array}{cc}\mathbf{\gamma }& 0\\ -\mathbf{I}& 0\end{array}\right),$$

$$\mathbf{X}(t^{\prime })=e^{-\hat{\mathbf{\gamma }}(t^{\prime }-t)}\mathbf{X}(t)+{\int }_t^{t^{\prime }}\text{d}s{e}^{\hat{\mathbf{\gamma }}(s-t^{\prime })}\left(\mathbf{F}(\mathbf{X}(s))+\mathbf{E}(s)\right).$$

$$\hat{\mathbf{\gamma }}=\frac{\mathbf{\gamma }}{2}\otimes \mathbf{I}-\frac{\mathbf{I}}{2}\otimes {\mathbf{\sigma }}_x+i\frac{\mathbf{I}}{2}\otimes {\mathbf{\sigma }}_y+\frac{\mathbf{\gamma }}{2}\otimes {\mathbf{\sigma }}_z.$$

$${\mathbf{p}}_{n+1}$$

$${\mathbf{q}}_{n+1}$$

$$\mathbf{\Gamma }$$

$$\mathbf{\Theta }$$

$$\tilde{\mathbf{\eta }}$$

$$\widetilde{\stackrel{~}{\mathbf{\eta }}}$$

$$Z\simeq \frac{1}{(2\pi \hslash {)}^f}\int {\text{d}}^f\mathbf{p}{\text{d}}^f\mathbf{q}{e}^{-\tau H_M(\mathbf{p},\mathbf{q})},$$

$$H_M(\mathbf{p},\mathbf{q})$$

$$\mathcal{L}$$

$$\mathbf{\gamma }={\mathbf{\gamma }}^{\text{BO}}+{\mathbf{\gamma }}^{\text{harm}}.$$

$${\mathcal{L}}^{\text{harm}}$$

$${\mathcal{L}}^{\text{BO}}$$

$$e^{-\mathcal{L}\delta t}\simeq e^{-{\mathcal{L}}^{\text{BO}}\delta t/2}{e}^{-{\mathcal{L}}^{\text{harm}}\delta t}{e}^{-{\mathcal{L}}^{\text{BO}}\delta t/2}.$$

$$\hat{\mathbf{\gamma }}=\left(\begin{array}{cc}{\mathbf{\gamma }}^{\text{harm}}& \mathbf{K}\\ -\mathbf{I}& 0\end{array}\right),$$

$$K_{ih}^{(j)(k)}={\tilde{\omega }}_M^2{\delta }_{ih}\left(2{\delta }^{(j)(k)}-{\delta }^{(j)(k-1)}-{\delta }^{(j)(k+1)}\right).$$

$$\hat{\mathbf{\gamma }}=\frac{{\mathbf{\gamma }}^{\text{harm}}}{2}\otimes \mathbf{I}+\frac{\mathbf{K}-\mathbf{I}}{2}\otimes {\mathbf{\sigma }}_x+i\frac{\mathbf{K}+\mathbf{I}}{2}\otimes {\mathbf{\sigma }}_y+\frac{{\mathbf{\gamma }}^{\text{harm}}}{2}\otimes {\mathbf{\sigma }}_z.$$

$${\mathbf{p}}_{n+1}$$

$${\mathbf{q}}_{n+1}$$

$${\mathbf{E}}_{\text{int}}=\left(\begin{array}{c}\mathbf{\Gamma }\tilde{\mathbf{\eta }}\\ \mathbf{\Theta }\widetilde{\stackrel{~}{\mathbf{\eta }}}\end{array}\right)$$

$$⟨{\mathbf{E}}_{\text{int}}^T{\mathbf{E}}_{\text{int}}⟩=\underset{-\delta t}{\overset{0}{\int }}dt{e}^{\hat{\mathbf{\gamma }}t}\hat{\mathbf{\alpha }}{e}^{{\hat{\mathbf{\gamma }}}^{†}t},$$

$$\hat{\mathbf{\alpha }}=\left(\begin{array}{cc}\mathbf{\alpha }& 0\\ 0& 0\end{array}\right).$$

Full text

I Introduction

In spite of about half a century of numerical investigations aimed at providing a quantitative and reliable microscopic description of liquid water, there is still a difficulty in describing accurately proton hopping inside the strong and well structured hydrogen bond network giving water its numerous peculiar features. Since water is the solvent of most chemical and biological processes occurring on Earth, one of the long-term goals is to explicitly take into account, at a high accuracy level, hydrogen-bond driven phenomena in the treatment of species interacting with this solvent. The solvation and transport properties of hydronium ($\text{H}_{3}\text{O}^{+}$) and hydroxide (HO-) ions in aqueous environment have enormous impacts on areas ranging from aqueous acid-base chemistryEigen (1964); Mohammed et al. (2005, 2007); Hassanali et al. (2011); Lee and Rasaiah (2011); Agmon et al. (2016), enzymatic proton transferRiccardi et al. (2006, 2010), as well as proton transfer in biological channelsWraight (2006), through fuel cell membranesDecoursey (2003); Freier et al. (2011) and on ice surfaces facilitating atmospheric reactionsGertner and Hynes (1998). Recent progress in our understanding of proton transport in water has been reviewed in Refs. Marx, 2006 and Hassanali et al., 2014.

Due to a variety of interactions present in water such as weak and dispersive dipolar interactions that compete during the formation and the breaking of hydrogen bonds inside a strongly polar network, an accurate description of the Potential Energy Surface (PES) is required to fully understand the microscopic properties of water. From a theoretical point of view, liquid water has been studied through Molecular Dynamics (MD) techniques such as the empirical force fields methodsJorgensen et al. (1983); Sprik (1991); Mahoney and Jorgensen (2000) and ab initio molecular dynamics (AIMD)Silvestrelli and Parrinello (1999); Grossman et al. (2004); VandeVondele et al. (2005); Jonchiere et al. (2011); Carloni et al. (2002); Friesner and Dunietz (2001). Because of its good scalability with the system size which is suited to study liquids, Density Functional Theory (DFT) has been the most employed technique to perform AIMD simulations of bulk water. Nevertheless, DFT-based calculations performed with Generalized Gradient Approximation (GGA) functionals give rise to an overstructuration of water: melting point temperatureYoo et al. (2009) and oxygen-oxygen radial distribution functionLee and Tuckerman (2007) are for instance poorly reproduced. Some improvements by including Van der Waals corrections such as the BLYP-DCACP, BLYP-DG2 or PBE0+TS-vdW(SC) functionals have been tested and proven to get results closer to experimentsLin et al. (2012); Santra et al. (2015) but these approaches suffer from a lack of generality and transferability. Gillan et al. recently established a numerical scoring system to evaluate the quality of different exchange-correlation functionals for a wide range of water systemsGillan et al. (2016).

As an alternative to DFT-AIMD, very accurate calculations using embedded-fragment second-order Møller-Plesset perturbation theory (MP2) have been carried out and gave satisfactory results despite the shortness of the MD trajectories (6 ps of equilibration and 5 ps of production run)Willow et al. (2015). Perturbative theories such as MP2/MP4 or the most advanced quantum chemistry techniques such as Coupled Cluster (CC) calculations are indeed suffering from a poor scalability with the number of particles $N$ in the system.

To overcome this major limitation, many-body representations of water and other aqueous systems appear to be a possible route to describe accurately proton transfer properties at a cost marginally higher than that associated with empirical force fields. The many-body approach has witnessed tremendous recent advances since numerous polarizable potentials have been proven to be closely accurate to CCSD(T) reference.Cisneros et al. (2016) Some of them, namely the MB-pol Babin et al. (2013, 2014); Medders et al. (2014), provide a correct description of water from the gas to condensed phase, outperforming all models currently available. Building upon this formalism, several many-body representations have also been introduced for protonated water.Pinski and Csányi (2014); Yu and Bowman (2016)

While keeping an explicit description of the electronic structure, Quantum Monte Carlo (QMC) techniques constitute a very interesting alternative to study aqueous systems, since the scalability of the method ($N^{3-4}$) makes the simulations of very large systems computationally affordable with a much greater accuracy than DFT in most cases. Recently, it has been demonstrated that QMC techniques provide results as accurate as the basis set converged CCSD(T) for small neutral or charged water clustersGillan et al. (2012); Benedek et al. (2006). Consequently, QMC can now be used as a benchmark method and certainly benefits from its intrinsically parallel formulation in modern supercomputers. The price to pay here is the systematic statistical incertitude arising from the stochastic nature of the QMC approach, which can however be reduced by generating longer Markov chains as the error is inversely proportional to the square root of the number of QMC generations.

Some strategies have thus been developed to incorporate the intrinsic QMC noise into a classical statistical mechanics framework of nuclei at finite temperature $T$. For instance, the Coupled Electron Ion Monte Carlo (CEIMC) method relies on the Born-Oppenheimer (BO) approximation for treating finite temperature ions coupled with ground state electrons. The Boltzmann distribution function of the ionic degrees of freedom is sampled at fixed temperature via a Metropolis MC simulation based on the electronic energies computed during independent ground state QMC calculationsPierleoni et al. (2004); Pierleoni and Ceperley (2006). Another strategy is to resort to a Langevin MD approach to correlate the noise driving the Langevin Dynamics (LD) of the system by the QMC forces covariance. In that case, the wave function parameters are optimized by Variational Monte Carlo (VMC) calculations along the MD path such that the electronic solution has always energies and forces as close as possible to the true BO surfaceAttaccalite and Sorella (2008); Mazzola et al. (2014).

Last year, Zen and coworkers applied the latter method to perform the very first QMC-based MD simulation of bulk liquid water at ambient conditions using 32 water molecules in a cubic cell with periodic boundary conditionsZen et al. (2015). By obtaining a better position of the first peak of the oxygen-oxygen radial distribution function compared with most advanced functionals, they demonstrated the ability of QMC to tackle this kind of problems. However, the authors had to make some approximations: they expanded the wave function over a small basis set and, more importantly, a classical description of the nuclei was adopted.

Indeed, it is well known that Nuclear Quantum Effects (NQE) play a crucial role in the description of water or ice by deeply affecting the radial distribution functions and distorting hydrogen bonds because of quantum disorderMorrone and Car (2008); Vega et al. (2010); Pamuk et al. (2012); Kong et al. (2012); Giberti et al. (2014); Lobaugh and Voth (1996); Markovitch et al. (2008). In this paper, we propose to extend Zen’s pioneering work by including a nuclear quantum description within the QMC-driven dynamics. This is achieved within a Path Integral Langevin Dynamics (PILD) approach which, to the best of our knowledge, has never been used in the case of non-deterministic forces.

On the one hand, we derive an original algorithm to integrate the Langevin equations of motion for quantum particles. The designed algorithm is very general since it can be used to propagate both deterministic and stochastic forces. On the other hand, within our framework we show that the inclusion of the quantum nature of the nuclei in the QMC-driven MD does not slow down the simulation with respect to the classical counterpart. This is thanks to an effective average of wave function parameters in the polymer representation of quantum nuclei, which relies upon a bead-grouping technique. Therefore, the cost of the wave function optimization is not increased by the inclusion of the additional quantum degrees of freedom.

To test our methodological development, we apply it to $\text{H}_{5}\text{O}_{2}^{+}$, namely the Zundel ionZundel and Metzger (1968), widely used as benchmark system. Indeed, it is the simplest system to exhibit a non trivial proton transfer and its reduced size makes a comprehensive and systematic study of the problem easier. More importantly, there is a huge amount of dataWei and Salahub (1994); Tuckerman et al. (1997); Vuilleumier and Borgis (1998a, b, 2000); Valeev and Schaefer (1998); Auer et al. (2000); Wales (1999); Xie et al. (1994); Headrick et al. (2005); Park et al. (2007) to compare our results with, because the description of excess proton in water has been massively studied both theoretically and experimentally in the last fifty years. We can cite for instance the extremely accurate Potential Energy Surface (PES) generated by Bowman and coworkers from almost exact CCSD(T) calculationsHuang et al. (2005), the MS-EVB methodVuilleumier and Borgis (1998a, b); Schmitt and Voth (1998); Day et al. (2002), or more recently the LEWIS model developed by Herzfeld.Kale and Herzfeld (2012) Moreover, because of the great importance of the NQE in the Zundel cation, the latter is almost always used to test and validate new approachesKapil et al. (2016).

The paper is organized as follows. In Section II, we introduce the theoretical framework by first reviewing the algorithms for classical Langevin simulations with QMC forces (Sec. II.1). In Sec. II.2 we extend the classical description of nuclei by introducing an original and efficient algorithm to integrate the equations of motion in the nuclear quantum case with or without noisy forces. A study of the numerical stability of the new path integral integrator is done in Sec. II.3 for deterministic forces, where a direct comparison with existing PILD algorithms is possible. The inclusion of QMC forces is discussed in Sec. III. The general form of the QMC wave function for the electronic part is described in Sec. III.1, the bead-grouping approximation is introduced in Sec. III.2, while the QMC-noise correction scheme is detailed in Sec. III.3, which yields an unbiased sampling of the canonical quantum partition function even in the presence of stochastic errors. Sec. III.4 reports on the numerical stability analysis of the new integrators with noisy QMC forces. In Section IV, we show the results obtained for the Zundel ion at low and room temperatures by taking into account, or not, the NQE and demonstrate the excellent agreement between the noisy QMC PES and the very accurate deterministic CCSD(T) one, which validates our method. Finally, in Section V, we provide some concluding remarks and we discuss the possible future applications of the algorithms developed in this paper.

II Methods

In this Section, we derive the methodology to sample the canonical partition function corresponding to Hamiltonians including either classical or quantum nuclei. This is achieved using the generalized Langevin dynamics, which has been demonstrated to be compatible with the path integral approach to include a quantum description of the nucleiCeriotti et al. (2011). Thanks to the stochastic nature of this category of non-Newtonian dynamics, the Langevin formulation of the equations of motion can be extended to the QMC case, where the forces are intrinsically noisyAttaccalite and Sorella (2008); Mazzola et al. (2014).

II.1 Classical Langevin Dynamics: Classical Momentum-Position Correlator (CMPC) algorithm

Let us consider a system constituted by $N$ classical particles described by the Hamiltonian

[TABLE]

where $m_{i}$ are the physical masses, $p_{i}$ and $q_{i}$ are the atomic momenta and positions, respectively, and $V$ is the PES in which the system evolves. Hereafter, we are going to use transformed variables, more convenient to manipulate, achieved by the following mass scaling:

[TABLE]

where $\left\{q_{i}^{0},p_{i}^{0},\eta_{i}^{0},f_{i}^{0}\right\}_{i=1,..,3N}$ are the original coordinates and $\left\{q_{i},p_{i},\eta_{i},f_{i}\right\}_{i=1,..,3N}$ denote the transformed ones. After applying the above transformation, the corresponding second order Langevin dynamics reads:

[TABLE]

where $\mathbf{f}=-\bm{\nabla}_{\mathbf{q}}V(\mathbf{q})$ is the applied force deriving from the PES, and the vectors are $3N$-dimensional. To ensure that the equilibrium distribution generated by this dynamics is canonical, the r.h.s of Eq. (3) includes a $3N\times 3N$ friction matrix $\bm{\gamma}$ and a random force vector $\bm{\eta}$. The latter controls the temperature $T$ of the system via the Fluctuation-Dissipation Theorem (FDT)Kubo (1966)

[TABLE]

with $k_{B}$ the Boltzmann constant, and

[TABLE]

is the force covariance matrix which simply reduces to a diagonal form for white noiseFixman (1986); Izaguirre et al. (2001). The dynamics generated by the stochastic differential Eqs. (3) and (4) is Markovian, as the noise fluctuations are locally correlated in time. In the case the ionic forces are computed by QMC methods (as we will see in Sec. III.3), the covariance matrix has non-trivial off-diagonal elements, and $\bm{\eta}$ fluctuations are correlated by the QMC statistics. Thus, the Langevin approach can naturally deal with noisy QMC forces, as already shown by Sorella and coworkers in the classical MD frameworkAttaccalite and Sorella (2008); Mazzola et al. (2014); Zen et al. (2015).

The basic integration scheme developed by Attacalite *et al.*Attaccalite and Sorella (2008) to perform QMC-driven MD simulations at finite temperature has been further developed in Refs. Mazzola et al., 2014 and Zen et al., 2015. In the latter scheme, dubbed as Classical Momentum-Position Correlator (CMPC) algorithm, the Langevin dynamics is driven by a correlated noise affecting both momenta and positions. To set up the formalism that we will be exploiting in the quantum path integral case, we provide here an alternative derivation, based on the joint momentum-position coordinates

[TABLE]

where $\mathbf{X}$ is a $6N$-dimensional vector. Analogously, we extend the definition of the ionic and random force vectors to be $\mathbf{F}=\begin{pmatrix}\mathbf{f}\\ \mathbf{0}\end{pmatrix}$ and $\mathbf{E}=\begin{pmatrix}\bm{\eta}\\ \mathbf{0}\end{pmatrix}$, respectively. In this extended basis, Eqs. (3) and (4) can be rewritten into a generalized SDE, which reads as

[TABLE]

where in this notation the matrix $\hat{\bm{\gamma}}$ represents a generalized friction that couples both momenta and positions by:

[TABLE]

with the “physical” friction $\bm{\gamma}$ being the same $3N\times 3N$ matrix introduced in Eq. (3), and $\mathbf{I}$ is the identity matrix. The formal solution of Eq. (8) is provided by

[TABLE]

Using joint coordinates would be of little use, if we were not be able to evaluate the exponential $e^{-\bm{\hat{\gamma}}\delta t}$ in a closed analytic form. This is possible because the block matrix $\bm{\hat{\gamma}}$ can be more conveniently rewritten in terms of Pauli matrices $\bm{\sigma}_{x}$, $\bm{\sigma}_{y}$, $\bm{\sigma}_{z}$, as follows:

[TABLE]

Then, the exponentiation can be straightforwardly obtained by using standard Pauli matrices algebra, and the solution can be given in a closed form:

[TABLE]

where the time evolution has been discretized with time step $\delta t$, and the subscripts refer to the corresponding time slice, such that $\mathbf{p}_{n}=\mathbf{p}(t_{n})$, $\mathbf{q}_{n}=\mathbf{q}(t_{n})$, and $\mathbf{f}_{n}=\mathbf{f}(\mathbf{q}(t_{n}))$. In the above Equations, the other symbols are defined as

[TABLE]

We immediately notice that the CMPC algorithm differs from those previously developed in literatureSkeel and Izaguirre (2002); Ricci and Ciccotti (2003); Vanden-Eijnden and Ciccotti (2006); Bussi and Parrinello (2007); Melchionna (2007); Grønbech-Jensen and Farago (2013), since momenta and positions are propagated simultaneously in a single iteration thanks to the use of momentum-position correlation matrices. In particular, according to Eqs. (12) and (13), not only the momenta but also the positions are affected by the integrated Langevin noise, $\bm{\tilde{\eta}}$ and $\bm{\tilde{\tilde{\eta}}}$ whose properties are detailed in Appendix A.

It is interesting to note that the momentum-position correlator formalism correlates the integrated noise, even without dealing with non-diagonal $\bm{\alpha}$ matrices, as we will do in the QMC case (Sec. III.3). We also remark that in order to derive the integrated equations of motion (12) and (13), we disregarded the $\mathbf{q}$-time dependence of $\mathbf{f}$ in Eq. (10) in the time interval $\delta t$. For deterministic forces, where $\bm{\alpha}$ is position independent, it is the only approximation left in the discretized classical Langevin dynamics driven by Eqs. (12) and (13).

For deterministic $\mathbf{f}$, $\gamma_{ij}=\gamma_{\textrm{BO}}\delta_{ij}$, while the definition of $\bm{\gamma}$ is more general for noisy QMC forces (see Sec. III). Once $\bm{\gamma}$ and $\delta t$ have been set, the numerical evolution is performed according to Eqs. (12) and (13) in the frame which diagonalizes $\bm{\gamma}$.

II.2 Path Integral Ornstein-Uhlenbeck Dynamics (PIOUD)

As we mentioned in the Introduction, it is well known that NQE are extremely important, in liquid water and aqueous species such as the Zundel ion, because of the light mass of the hydrogen atom. Several approaches to deal with NQE have been proposed so far, such as the Quantum Thermal Bath (QTB)Dammak et al. (2009); Ceriotti et al. (2009a); Brieuc et al. (2016a, b), but we have chosen here to focus on the Path Integral (PI) framework for the following reasons. First, the formulation of the theory is intrinsically parallel and perfectly compatible with the QMC multi-walker implementation. Second, the PI relies on approximations which are systematically improvable and the exact result can be recovered with controlled extrapolations. Third, this approach can be easily combined with the LD leading to the PILD method, for which very efficient algorithms already existCeriotti et al. (2010a); Liu et al. (2016).

We first recall the basic formalism of the PIMD method (see Ref. Ceriotti et al., 2010a for further details), before introducing our approach to propagate the corresponding equations of motion. In PI theory, the quantum partition function $Z=\text{Tr}\left[e^{-\beta H}\right]$ (with $\beta=1/k_{B}T$) is usually approximated after a standard Trotter factorizationTrotter (1959) by

[TABLE]

with $f=MN$ and $\tau=\beta/M$. We have replaced here the true quantum particles by fictitious classical ring polymers consisting of $M$ replicas (beads) of the system connected to each other by harmonic springs evolving in the imaginary time $\tau$. Without loss of generality, we can extend the definition of classical vectors to the quantum case, by including all $M$ replicas in the string, such that $\mathbf{q}\equiv\left(\mathbf{q}^{(1)},\ldots,\mathbf{q}^{(i)},\ldots,\mathbf{q}^{(M)}\right)^{T}$ and it becomes $3NM$-dimensional. Bold vectors will automatically include all replicas in the quantum case. If interpreted classically, the partition function $Z$ in Eq. (18) describes a system at the effective temperature $MT$. The Hamiltonian $H_{M}$ corresponding to this quantum-to-classical isomorphism reads as

[TABLE]

written in mass-scaled coordinates (2), with $\tilde{\omega}_{M}=M/\beta\hbar$, and subjected to the ring boundary condition $q_{i}^{(0)}\equiv q_{i}^{(M)}$. Thus, the PIMD method is a way to calculate quantum mechanical properties using the recipes from classical statistical mechanics with a modified Hamiltonian containing an additional quantum kinetic term.

Various choices are available for the thermostat, ranging from the more conventional Nosé-Hoover chainNosé (1984); Hoover (1985) to stochastic thermostatting in both its simple path integral Langevin equation (PILE) formCeriotti et al. (2010a), which fulfills the FDT, and its Generalized Langevin Equation (GLE) flavorsCeriotti et al. (2010a, b, 2011), where the noise is colored to simulate quantum effects and accelerate the convergence of the PIMD with the number of beads, at the price of breaking the validity of the FDT. Here, we keep fulfilling the FDT. However, instead of using the PILE integration algorithmCeriotti et al. (2010a), we will extend the CMPC algorithm, described in Sec. II.1, to the quantum case. As already mentioned, this framework will allow us to incorporate noisy QMC forces in the equations of motion, by using an appropriately tailored Langevin noise without breaking the FDT, as it is explained in Sec. III. Moreover, in our quantum algorithm, which makes use of a Trotter breakupTrotter (1959) between the harmonic and the physical modes, the quantum harmonic part is integrated exactly together with the Langevin thermostat for the harmonic frequencies. Indeed, the corresponding Langevin equations of motion are linear. It is then possible to exploit the mapping onto an Ornstein-Uhlenbeck process (OUP), which is exactly integrable.

To derive our algorithm, we employ the Liouvillian approach, that is used to integrate the Fokker-Planck differential equation, describing the evolution of the ensemble probability density $\rho(\mathbf{p},\mathbf{q},t)$Bussi and Parrinello (2007). Since the Fokker-Planck is a linear differential equation, this quantity evolves formally as $\rho(t+\delta t)=e^{-{\cal L}\delta t}\rho(t)$, with ${\cal L}$ being the Liouville operator, defined as

[TABLE]

in mass-scaled coordinates, with $\partial_{q_{i}}\equiv\frac{\partial}{\partial_{q_{i}}}$. ${\cal L}$ is built upon the Hamiltonian propagation, driven by the first two terms, and the Langevin thermostat, represented by the last two, deriving from the Fokker-Planck equation. For the sake of readability, $i$ and $j$ run here over all particles and beads indexed together. In Eq. (20), ${\cal F}_{i}\equiv f^{\textrm{BO}}_{i}+f^{\textrm{harm}}_{i}$ is the generalized force, comprising the BO and harmonic contributions, where $f^{\textrm{BO}}_{i}\equiv-\partial_{q_{i}}{\tilde{V}}$ and ${\tilde{V}}=\sum\limits_{j=1}^{M}V(q_{1}^{(j)},..,q_{3N}^{(j)})$.

The quantum-to-classical isomorphism Hamiltonian in Eq. (19) includes very different energy scales. To cope with this, we would like to split the Liouvillian in Eq. (20) into just two operators, one containing only the physical (BO) modes, the other depending exclusively on the fictitious (harmonic) modes. To do so, we first separate the friction matrix into two contributions:

[TABLE]

We can then rewrite the total Liouvillian as the sum of two terms, ${\cal L}={\cal L}^{\textrm{BO}}+{\cal L}^{\textrm{harm}}$, where

[TABLE]

in such a way that we can break up the evolution operator via Trotter factorizationTrotter (1959). In order to preserve the reversibility of the process and reduce at most the time step error, we approximate the propagation using a symmetric formTuckerman et al. (1992); Sexton and Weingarten (1992):

[TABLE]

The equations of motion corresponding to ${\cal L}^{\textrm{harm}}$ are those in Eq. (8), provided $\mathbf{X}$ is now interpreted as a $6NM$-dimensional vector, the thermal noise $\bm{\eta}$ lives in the $3NM$-dimensional space, and the generalized $\bm{\hat{\gamma}}$ must be redefined in order to include the harmonic couplings between the beads. $\bm{\hat{\gamma}}$ is now a $6NM\times 6NM$ matrix, which reads as

[TABLE]

where $\mathbf{K}$ is $3NM\times 3NM$ matrix defined as follows:

[TABLE]

In the above definition, we have used lower indices to indicate the particle components, the upper ones (in parenthesis) indicate the bead components, while $\delta_{ij}$ ($\delta^{(i)(j)}$ for the bead indices) is the usual Kronecker delta. The $\mathbf{K}$ matrix is diagonal in the particle sector, as the harmonic springs in the fictitious $H_{M}$ of Eq. (19) couple different replicas only for the same particle components. Cyclic conditions are implicitly applied in Eq. (26) to the matrix boundaries in the bead sector (i.e. $(0)=(M)$), as the polymers are necklaces.

The most notable difference of the ${\cal L}^{\textrm{harm}}$ equations of motion with respect to those of Eq. (8) is that $\mathbf{F}=\mathbf{0}$, as the BO component of the total force has been loaded in the ${\cal L}^{\textrm{BO}}$ Trotter factors. Therefore, the ${\cal L}^{\textrm{harm}}$ equations are linear in both $\mathbf{p}$ and $\mathbf{q}$, thus exactly solvable in an analytic closed form. They describe an OUP, and their algebra and properties have been exploited in the GLE frameworkCeriotti et al. (2010b), where they are used to propagate external auxiliary variables together with the physical momenta in order to generate a colored noise with a corresponding effective memory kernel. In the present algorithm we are going to use the OUP to integrate directly the SDE for both $\mathbf{p}$ and $\mathbf{q}$, together with the Langevin noise, namely without further splitting the Langevin thermostat in ${\cal L}^{\textrm{harm}}$.

As we have done in Eq. (11) for the classical case, we expand $\bm{\hat{\gamma}}$ of Eq. (25) in Pauli matrices, as follows:

[TABLE]

Then, Eq. (27) and $\left[\mathbf{K},\bm{\gamma}^{\textrm{harm}}\right]=0$ (condition fulfilled, being $\bm{\gamma}^{\textrm{harm}}$ bead independent) allow one to evaluate $e^{\bm{\hat{\gamma}}\delta t}$ in a closed analytic form for each (upper and lower) block component of the SDE formal solution in Eq. (10). The related algebra is quite tedious and we refer the reader to the Appendix C. The resulting integrated equations of motions lead to the following Markov chain:

[TABLE]

with the integrated $6NM$-dimensional noise ($\mathbf{E}_{\textrm{int}}$) which is

[TABLE]

The expressions for $\bm{\Lambda}$, $\bm{\Gamma}$ and $\bm{\Theta}$ matrices are quite complex, so they are given in Appendix C. Similarly to the classical case, we can compute the noise correlation matrix by using its definition:

[TABLE]

where

[TABLE]

Once again, we refer the reader to the Appendix C for the analytic expression of the noise correlators in the quantum case.

The propagator $e^{-{\cal L}^{\textrm{harm}}\delta t}$ is applied exactly, and the resulting evolution is given by the coupled equations (28) and (29) solved in the normal modes representation. $\bm{\gamma}^{\textrm{harm}}$ is defined to be the optimal damping for the dynamics driven by a harmonic $H$ in its normal modes representationTassone et al. (1994); Ceriotti et al. (2010a), namely

[TABLE]

for each harmonic eigenmode $k=1,\ldots,M$, where

[TABLE]

In other words, the thermostat controlled by the friction $\gamma_{0}$ is mainly applied to the centroid rototranslational modes which would not be thermalized otherwise because their frequency is zero. Contrary to the other free quantum harmonic modes, the optimal damping $\gamma_{0}$ is not general and has to be optimized for the system under study. A reasonable value of this parameter, also present in the PILE method, can be found by performing short preliminary simulations with any accurate force field Góra et al. (2014); Babin et al. (2012, 2013, 2014).

The equation of motion corresponding to ${\cal L}^{\textrm{BO}}$ is the one in Eq. (12) of Ref.Luo et al., 2014. Indeed, only the momenta are evolved in ${\cal L}^{\textrm{BO}}$. Therefore, no MPC algorithm is needed, and the resulting equations are equal to those of the simple classical Langevin algorithm introduced in Ref.Attaccalite and Sorella, 2008, restricted to $\mathbf{p}$, which in the present case is a $3NM$-dimensional vector. The corresponding $\bm{\gamma}^{\textrm{BO}}$ is taken block diagonal according to the bead index, and possibly bead dependent. This is very useful for QMC noisy forces, as the statistics of the QMC force covariance matrix is genuinely bead dependent, because the wave function is sampled independently for each bead (see Sec. III). In the QMC case, $\bm{\gamma}^{\textrm{BO}}$ is fundamental to correct the intrinsic QMC noise affecting the BO forces. However, for deterministic BO forces, the corresponding $\bm{\gamma}^{\textrm{BO}}$ can be taken equal to zero, and the corresponding time-discretized solution will result in a simple velocity update $\mathbf{p}_{n+1}=\mathbf{p}_{n}+\frac{\delta t}{2}\mathbf{f}_{n}$. In fact, any additional damping in the BO modes will turn into a slower algorithmic diffusion, as the thermalization is already guaranteed by the Langevin thermostat integrated in ${\cal L}^{\textrm{harm}}$. In this limit, ${\cal L}^{\textrm{BO}}$ will reduce to ${\cal L}_{\textrm{V}}$, using the same notations as in Ref. Ceriotti et al., 2010a.

For deterministic forces, Ceriotti and coworkersCeriotti et al. (2010a) developed a different algorithm (dubbed PILE), where the main difference with respect to ours is that ${\cal L}^{\textrm{harm}}$ is split into ${\cal L}_{0}$, which propagates the motion of the harmonic oscillators in their normal modes representation, and ${\cal L}_{\gamma}$ (the part proportional to $\gamma_{ij}^{\rm harm}$ in Eq. (22), ${\cal L}_{0}$ being the remaining part in the first line of the same equation), which corresponds to the Langevin thermostat acting on the normal modes. The total Liouvillian propagator proposed in Ref. Ceriotti et al., 2010a is:

[TABLE]

Apart from the additional Trotter breakup, in the above sequence of operators the Langevin thermostat on the normal modes is the outermost part of the Liouvillian decomposition, for a better control of the target temperature. Our algorithm in Eq. (24) performs better than the solution proposed in Eq. (35), since in our case ${\cal L}^{\textrm{harm}}={\cal L}_{0}+{\cal L}_{\gamma}$ corresponds to an exact propagator including both the coordinates change and the thermalization of the free quantum ring polymers. As a consequence, one Trotter factorization is saved which in principle enables us to work with a larger time step $\delta t$ or more quantum replicas $M$. Indeed, the Trotter break-up is obviously exact if the two operators involved commute. The split of $\mathcal{L}^{\textrm{harm}}$ proposed by Ceriotti is a very unfortunate choice as the commutator of $[\mathcal{L}_{0},\mathcal{L}_{\gamma}]$ is dominated by the term $\gamma_{ij}^{\textrm{harm}}f^{\textrm{harm}}_{j}\partial_{p_{i}}$ which diverges as $M^{3}$ for large $M$, as it can be easily seen by the exact expression of the harmonic forces and the choice of the optimal friction. This argument suggests that our method, involving commutators at most diverging as $M^{2}$, should have a time step error much better behaved for large $M$.

These algorithmic improvements will be tested and discussed in the next section.

II.3 Algorithm stability with deterministic forces

We test the robustness of the PIOUD algorithm detailed in Sec. II.2 by comparing it with PILE. We perform PILD simulations on the Zundel ion with almost exact deterministic forces which are simply computed by finite differences of the CCSD(T) PES provided by Huang and coworkersHuang et al. (2005). Since our only aim here is to test the different integration schemes, the number of MD iterations $N_{\textrm{iter}}=100000$ is quite small to visit the entire phase space but sufficient to identify some possible weakness of the propagators.

A robust PILD algorithm must be stable with both large time step $\delta t$ and large number of quantum replicas $M$. Indeed, the collective modes of a large ring polymer become very stiff, so much more difficult to control, especially if the integration time step is large. Close values of the virial $T_{M,\textrm{vir}}(\mathbf{q})=\frac{N}{2\beta}+\frac{1}{2M}\sum\limits_{i=1}^{3N}\sum\limits_{j=1}^{M}\left(q_{i}^{(j)}-\bar{q_{i}}\right)\partial_{q_{i}^{(j)}}{\tilde{V}}$ and the primitive $T_{M,\textrm{pri}}(\mathbf{q})=\frac{3NM}{2\beta}-\frac{1}{M}\sum\limits_{j=1}^{M}\frac{1}{2}\left(\frac{q_{i}^{(j)}-q_{i}^{(j-1)}}{\hbar\tau}\right)^{2}$ kinetic energy estimators point that we are sampling properly both positions and momenta of all particles.

II.3.1 Stability with respect to the number of beads $M$

We report the results obtained for the PILE and PIOUD integration schemes at low temperature ($T=50$ K) and room temperature ($T=300$ K) in Figure 1 as a function of the number of quantum replicas $M$. In each case, we compare the key observables $T_{\textrm{vir}}$ and $T_{\textrm{pri}}$ (other useful quantities, such as the average temperature $T$, the algorithmic diffusion constant $D$, and the potential autocorrelation time $\tau_{V}$, are discussed in the Supplementary Information (SI)). We perform simulations working with an increasing number of beads $M=4,8,16,32,64,128$ and $256$, using a reasonable but not-so-small time step $\delta t=0.5$ fs.

The quantum kinetic energy estimators are known to be useful tools to determine the number of beads required at a given temperature to capture NQEs. In the case of the Zundel ion, $M=128$ quantum replicas are enough to fully recover the quantum kinetic energy at low temperature, whereas only $M=32$ beads are required at room temperature, according to Figure 1. We note that these values are consistent with those frequently used in the literature for this systemMarx et al. (1999); Spura et al. (2015); Suzuki et al. (2013).

The computational efficiency and the performances of the PIOUD propagator of Eq. (24) and the PILE algorithm are very close to each other, regarding the average simulation temperature or the algorithmic diffusion constant (see SI for more details). This is expected because the same normal mode transformation is applied in both approaches.

Nevertheless, when looking at the Figure 1, we notice that the PIOUD algorithm displays an almost perfect control of the kinetic energy operators at each temperature. On the contrary, the PILE algorithm exhibits a significant instability of the primitive energy for a large number of beads at both low and room temperature. The robustness of our approach mainly relies on a simultaneous control of both positions and velocities during the dynamics via momentum-position correlation matrices.

In summary, the PIOUD propagator appears to be more efficient than PILE as the positions are better controlled, which enables us to work with a larger number of quantum replicas at a fixed time step $\delta t$.

II.3.2 Stability with respect to the time step $\delta t$

Thanks to the preliminary analysis discussed in Sec. II.3.1, we know how many quantum replicas should be used to generate an efficient and converged PILD simulation of the Zundel ion at a given temperature. Moreover, stability for rather large time steps is crucial in the perspective of performing PILD calculations on larger systems using a PES evaluated by accurate but computationally demanding ab initio methods. In this Section, we report CC-PILD simulations at room temperature with $M=32$ beads using different values of $\delta t=0.1,0.2,0.3,0.5,0.75$ and $1$ fs. The behavior of the observables used to evaluate the algorithm efficiency are plotted in Figure 2 as a function of the time step at 300 K.

In order to quantify the bias induced by the time step error, we check the difference $|\langle T\rangle_{M,\textrm{vir}}-\langle T\rangle_{M,\textrm{pri}}|$ which gives us direct information on the accuracy of the positions sampling. Like in the previous tests, the PIOUD algorithm shows the smallest difference thanks to a good control of the primitive energy, due to a better integration scheme. The difference between these two kinetic energy estimators is quite spectacular at room temperature. Indeed, the fluctuation-dissipation contributions in ${\cal L}^{\textrm{harm}}$, related to damping and random forces in the dynamics, become more important as the temperature increases, while the BO forces are not so strongly affected by thermal effects.

The supplemental analysis based on $T$, $D$, and $\tau_{V}$ is reported in the SI. The algorithmic tests carried out here have been applied to a deterministic PES and exact analytic forces. We now wish to go further by applying the previous formalism to the case of stochastic PES and forces, such as the ones computed in a QMC framework, without losing computational efficiency.

III Extension to the stochastic case: correlating the noise by QMC

III.1 QMC-BO surface and electron wave function

After checking the performances of the integrators with deterministic forces, we shall show that the PIOUD algorithm can be easily extended to perform PILD simulations in a potential energy landscape $V(\mathbf{q})$ evaluated by VMC, namely:

[TABLE]

where $|\Psi_{\mathbf{q}}\rangle$ is the QMC wave function, that minimizes the expectation value of $H(\mathbf{q})$. In the QMC framework, the evaluation of $V(\mathbf{q})$ and its corresponding force estimators are intrinsically noisy. In Eq. (36) and in the following, we use the $\mathbf{q}$ subscript to make apparent that the VMC wave function has an explicit dependence on the nuclear coordinates. Using the BO approximation, the PES corresponds to the expectation value of the exact Hamiltonian applied to the wave function $|\Psi_{\mathbf{q}}\rangle$ written here as a Jastrow Antisymmetrized Geminal Power (JAGP) productCasula et al. (2004)

[TABLE]

The set $\left\{\mathbf{x}_{i}=(\mathbf{r}_{i},\sigma_{i})\right\}_{i=1,\dots,N_{\mathrm{el}}}$ represents spatial and spin coordinates of the $N_{\mathrm{el}}$ electrons.

On the one hand, the bosonic Jastrow factor $J_{\mathbf{q}}$ is a function of electron-electron separation and mainly describes dynamical correlations due to charge fluctuations and Van der Waals effects. It is often written as a product of one-body, two-body and three/four-body terms $J_{\mathbf{q}}=J_{1,\mathbf{q}}J_{2,\mathbf{q}}J_{3,\mathbf{q}}$. The one body term reads

[TABLE]

with $u(|\mathbf{r}-\mathbf{q}|)=\frac{1-e^{-b|\mathbf{r}-\mathbf{q}|}}{2b}$ and $b$ is a variational parameter. This form is chosen to satisfy the electron-ion Kato cusp conditions, used to deal appropriately with the diverging electron-nucleus Coulomb potentials ad short distances. In the Zundel ion Hamiltonian, our benchmark system, we keep the bare Coulomb potential only for hydrogen, while for the oxygen atoms we replace it with the Burkatzki-Filippi-Dolg (BFD) pseudopotentialBurkatzki et al. (2007), which is smooth at the electron-ion coalescence points. Thus, $J_{1,\mathbf{q}}$ is applied only to the hydrogen atom. Similarly, the electron-electron cusp conditions are dealt with by the two-body Jastrow factor

[TABLE]

where $u$ is a function of the same form as in Eq. (38), but with a different variational parameter. Finally, many-body correlations are included in the remaining part of the Jastrow factor

[TABLE]

with

[TABLE]

where $N_{\mathrm{basis}}$ is the number of $\Psi_{a,\mu}^{J}$, the Gaussian Type Orbitals (GTOs) constituting the primitive basis of each atom $a$. The Jastrow functional form has been developed on a primitive Gaussian basis set of O(3s,2p,1d) H(2s,1p) corresponding to 213 Jastrow parameters.

On the other hand, since our system is spin unpolarized, the Fermionic part of the wave function is expressed as an antisymmetrized product of the geminals or pairing (AGP) functions $\Phi_{\mathbf{q}}(\mathbf{x}_{i},\mathbf{x}_{j})$:

[TABLE]

The geminals are antisymmetric functions of two electron coordinates written as a product of a spatial symmetric part and a spin singlet

[TABLE]

where the $\phi_{\mathbf{q}}(\mathbf{r}_{i},\mathbf{r}_{j})$ can be expanded either over the primitive basis (O(5s5p2d) and H(4s2p) for the determinantal part) or the computationally more efficient hybrid orbitals $\bar{\Psi}^{AGP}$, also developed on GTOs, such that

[TABLE]

Indeed, one can apply here the geminal embedding schemeSorella et al. (2015) to obtain $N_{\mathrm{hyb}}<N_{\mathrm{basis}}$ geminal embedding orbitals (GEOs), to reduce the total number of parameters, that is almost proportional to the computational cost spent for the optimization of the wave function. Here, we use the 8O 2H hybrid basis set which has been proven to be the best compromise between accuracy and computational cost, due to a limited number (571) parameters describing the determinantal part of the wave functionDagrada et al. (2014). This point is very important, because in the present optimization methodsSorella (2001); Casula et al. (2004); Umrigar et al. (2007), the matrices involved in the iterative procedure have linear dimension equal to the number $p$ of parameters involved. The number of QMC sampling used to characterize stochastically a $p\times p$ matrix should be much larger than $p$, otherwise the matrix is biased if not rank deficient, and the QMC optimization methods no longer work. Thus, it is clear that the reduction of the number of parameters is very much needed at present, because it is basically proportional to the computational cost of QMC optimization.

In this paper we use the optimal wave function devised by Dagrada et al. for the same system at zero temperature; further technical details concerning the variational wave function can be found in Ref. Dagrada et al., 2014.

During the dynamics, the GTO exponents in both the Jastrow and the AGP part of the wave function are kept frozen to make the simulation stable. At each new ionic configuration, the wave function must be reoptimized. This is done by means of the Stochastic Reconfiguration (SR)Sorella (2001); Casula et al. (2004), or the optimization method with Hessian acceleration (SRH)Umrigar et al. (2007). As the ionic positions are smoothly connected to those of the previous MD time step, also the wave function parameters will evolve continuously. After each nuclear iteration, the wave function is optimized by few (five in our simulations of the Zundel ion) SR or SRH steps on the electronic parameters to ensure that the system is close enough to the BO surface, before continuing the propagation of the ion dynamics. Due to the continuity of the nuclear trajectories, the number of SR steps is significantly smaller than the one required for a wave function optimization from scratch.

III.2 Bead-grouping approximation

Another advantage of our framework is represented by the local, ion-centered, basis set. We start from the observation that the most relevant dependence of the wave function on the ionic positions $\mathbf{q}$ comes directly (and explicitly) from the basis set, while the dependence of the electronic variational parameters is generally weaker. This is particularly true in quantum MD calculations, where each particle is represented by an $M$-bead necklace. In the first approximation, a necklace could share the same wave function parameters, while the full $\mathbf{q}$ dependence is provided by the basis set only, a dependence coming from the atomic centers -different in each bead- defining the basis set at each different (quantum) time slice. This approximation can be very effective to reduce the computational burden, because the main drawback encountered in PIMD and PILD calculations with respect to their classical counterparts is the factor-of-$M$ increase of variational parameters ($M\times p$ if no approximation is employed). Conversely, if the number of variational parameters is restricted to the same number as the classical simulation, no computational overhead is expected to work with $M>1$ in QMC. Indeed, the evaluation of the ionic forces can be done with a number of samples inversely proportional to $M$, as the temperature in each time slice is increased by $M$ and the Langevin equations require statistically less accurate -but nevertheless unbiased- forces, so that increasing $M$ becomes essentially cost-free in QMCPierleoni and Ceperley (2006).

When instead forces are computed with deterministic methods - e.g. DFT- the computational burden is necessarily proportional to $M$, and therefore several techniques have been recently developed to decrease the number of evaluations of the ionic forces by about an order of magnitude without missing any significant NQE. This is done by achieving smart interpolations and groupings of different possible pathsCheng et al. (2016); Poltavsky and Tkatchenko (2016) or by applying a Ring Polymer Contraction (RPC) schemeMarkland and Manolopoulos (2008a, b); Fanourgakis et al. (2009). Otherwise, the number of quantum replicas can be reduced by working with a GLE including a colored noise mimicking the quantum fluctuations of the nuclei Ceriotti et al. (2009b, a, 2010b, 2011).

Although these methods could be effectively incorporated into our QMC framework, we are forced to explore another approach, because, as we have discussed, in QMC the problem is just the large number of variational parameters, and not the large value of $M$. We introduce therefore a method which takes advantage of the explicit wave function representation of the electronic problem. Indeed, each bead at each iteration has its own optimal wave function $|\Psi^{(k)}_{\mathbf{q}}\rangle$, for $k=1,\dots,M$, which minimizes the variational energy at the nuclear configuration $\mathbf{q}^{(k)}$. Consequently, we need to find the best variational parameters set \bm{\lambda}^{(k)}=\bigl{\{}g_{\mu,\nu}^{a(k),b(k)},\lambda_{\mu,\nu}^{a^{(k)},b^{(k)}},b^{(k)}\bigr{.}, orbital contraction coefficients,\bigl{.}\ldots\bigr{\}} for each wave function. Despite the availability of efficient QMC optimization algorithms, that task is still computationally demanding and would limit the application of our method to very small systems. To overcome this major difficulty, we exploit the local nature of the Gaussian basis sets used in the expansion of both the Jastrow and AGP factors. As anticipated, we make the approximation of defining $N_{\mathrm{groups}}$ groups of neighboring beads and constraining the wave function parameters to be equal for all beads in the same group. Since a group shares the same parameters, the corresponding energy gradients are then averaged over the quantum replicas constituting the group. In this way, we improve the statistics by a factor of $M/N_{\mathrm{groups}}$. We obtain less noisy parameters even though the resulting wave function is not exactly optimized for each quantum replica. We mention that this is a controllable and systematically improvable approximation. Indeed, if one takes $N_{\mathrm{groups}}=M$, the electronic result is exact, whereas $N_{\mathrm{groups}}=1$ constitutes the roughest approximation. In the latter case, one performs a fully quantum dynamics with almost the same statistics as the one with classical nuclei. We are going to test and validate this approximation for the Zundel ion (see Sec. IV). In the case of the simple hydrogen molecule the approximation with $N_{\mathrm{groups}}=1$ allows one to recover $90\%$ of the Zero Point Energy (ZPE).

III.3 QMC ionic forces and noise correction

In our approach, the potential energy landscape and thus the ionic forces acting on each particles are evaluated by VMC. As already done Sec. II, we denote the $3NM$-dimensional force acting on all the $NM$ quantum replicas by

[TABLE]

where $\bm{\nabla}_{\mathbf{q}}$ is the gradient relative to the Cartesian coordinates $\mathbf{q}$ of all the images of the nuclei. Its expression is quite complex because it contains the implicit and explicit dependence of $V(\mathbf{q})$ on both the nuclear positions $\mathbf{q}$ and the electronic parameters set $\bm{\lambda}$. Computing these forces with finite variance and in a fast way is of paramount importance to make a QMC-based MD and PIMD possible. This has become feasible, thanks to recent improvements, such as the space warp coordinate transformation (SWCT)Filippi and Umrigar (2000), generalized to infinitesimal ion displacements via algorithmic differentiation (AD)Sorella and Capriotti (2010). Thanks to the AD, computing all components of the ionic forces is only four times more expensive than the cost of an energy calculation, which is computationally affordable.

Another issue related to the stochastic nature of this method is the control of the statistical noise introduced into the dynamics. The latter must be kept under control if we want to have an unbiased sampling of the phase space during the propagation of the trajectory. Fortunately, the methodology described in the previous sections is particularly suited for this situation, as it deals with the most general case, where the friction $\bm{\gamma}$ and covariance $\bm{\alpha}$ matrices are not diagonal. This corresponds to the inclusion of a correlated noise into the dynamics, in contrast to the more usual white noise case. By construction, we can now make the assumption that there exists also a QMC contribution to the covariance matrix which becomes position-time dependentAttaccalite and Sorella (2008)

[TABLE]

$\alpha_{\textrm{BO}}=2k_{B}T\gamma_{\textrm{BO}}$ is the white noise contribution and $\Delta_{0}$ is a tunable parameter which is mainly set to make the covariance (and thus the friction) matrix positive definite. $\bm{\alpha}$ and $\bm{\gamma}$ are linked by the FDT, while $\bm{\alpha}^{\textrm{QMC}}$ is the QMC-force covariance matrix defined by

[TABLE]

where $\delta f_{i}=f_{i}-\langle f_{i}\rangle$ is the fluctuation of the i-th ionic force. Previous works have established that the QMC covariance matrix is roughly proportional to the dynamical matrix (it is exactly proportional for harmonic forces). Therefore, it carries information on the vibrational properties of the systemLuo et al. (2014). As a consequence, one could take advantage of an apparent drawback of the QMC approach by using the intrinsic noise to drive a dynamics in the phase space with nearly optimal sampling.

However, to fulfill the FDT, the QMC noise has to be disentangled from the total Langevin noise in the following way:

[TABLE]

such that $\alpha^{\textrm{QMC}}_{ij}=\langle\eta_{i}^{\textrm{QMC}}\eta_{j}^{\textrm{QMC}}\rangle$, namely $\eta_{i}^{\textrm{QMC}}$ is the noise already present in the $i-$th QMC force component. Therefore, by using that the external noise is independent of the QMC one, we obtain

[TABLE]

The above equation is valid at the time instant $t$. However, in the schemes devised in Secs. II.1, and II.2, the equations of motion are discretized and integrated over a finite time step $\delta t$, and also the Langevin noise ${\tilde{\eta}}$ appearing in those equations is integrated over the same $\delta t$. This will slightly change the form in Eq. (49), leading to the following one:

[TABLE]

where ${\tilde{\alpha}}$ is the integrated noise correlator matrix, reported in Eq. (A) for the CMPC algorithm. $\bm{\alpha}^{\textrm{QMC}}$ does not need to be modified, as the expectation value of the variance-covariance QMC-force matrix in Eq. (47) is evaluated in this case during the time step $\delta t$ while sampling the electron coordinates. Since the r.h.s. of Eq. (50) is a positive definite matrix (thanks to the appropriately chosen $\Delta_{0}\geq\delta t$), it defines a corrected external noise which is compatible with the solution of the Langevin dynamics in both classical and quantum cases, as long as $|\mathbf{q}_{n+1}-\mathbf{q}_{n}|$ remains small, regardless of how large the friction is. As already mentioned, the proposed MPC schemes yield a correlated noise also affecting the conjugate variables (i. e. the nuclear positions $\mathbf{q}$). Therefore, the relation in Eq. (50) must be extended in the joint momentum-position coordinates. The full expression for the noise correction in the classical CPMC integration algorithm is detailed in Appendix B.

In the quantum PIOUD algorithm, the QMC noise correction must be applied exclusively in the BO step ${\cal L}^{\textrm{BO}}$, where only the momenta are evolved according to the QMC ionic forces. Therefore, the external noise acts only on the momentum sector, and its variance-covariance is simply given by Eq. 50. Contrary to the CPMC case, there is no need to extend this relation to the position sector. In the harmonic step ${\cal L}^{\textrm{harm}}$, the harmonic forces are noiseless, therefore no correction is needed, and the external noise coincides with the full Langevin noise ($\eta_{i}^{\textrm{QMC}}=0$). Its momentum and position components ($\bm{\tilde{\eta}}$ and $\bm{\tilde{\tilde{\eta}}}$, respectively) are correlated according to Eqs. 62 and 63, reported in Appendix C.

III.4 Algorithm stability with QMC forces

In the following, we shall present some calibration runs carried out with the PIOUD algorithm in the VMC framework, which will help us set the proper simulation parameters ($\gamma_{0}$, $\delta t$) and show the remarkable stability of the quantum dynamics even with noisy VMC forces for large $\delta t$ and large $M$. The results are obtained by performing QMC-PILD short tests (about $8.5$ ps of dynamics) of the Zundel ion in the gas phase at room temperature ($300$ K). Hereafter, the additional PIOUD parameter $\gamma_{\textrm{BO}}$ is set to zero, to avoid overdamping in the BO propagator. The Langevin thermostat in the ${\cal L}^{\textrm{BO}}$ part plays solely the role of correcting the BO dynamics for the intrinsic VMC noise. As we have already seen, for deterministic forces ${\cal L}^{\textrm{BO}}$ reduces to ${\cal L}^{V}$, i.e. it is a simple velocity step.

For each MD step we have evaluated forces and all the energy derivatives with $\approx 3.3\times 10^{5}$ MC samples, much larger than the total number of variational parameters ($p\approx 2.4\times 10^{3}$) of the wave function. This allows reaching an accuracy of $1.5$ mHa in the total energy per VMC energy minimization step. Between two MD steps, five QMC energy minimizations are performed with the Hessian (SRH) algorithm, in order to sample the correct BO surface. As already mentioned in Sec. III.1, all variational parameters are evolved during the dynamics, except for the GTO exponents, which are kept frozen.

The optimal input friction $\gamma_{0}$ is chosen to minimize the potential autocorrelation time $\tau_{V}=\frac{1}{\langle V^{2}\rangle-\langle V\rangle^{2}}\int\limits_{0}^{N_{\textrm{iter}}\delta t}dt\langle\delta V(0)\delta V(t)\rangle$ (with $\delta V=V-\langle V\rangle$ the fluctuation of the potential energy), and generate the most efficient phase space sampling during the dynamics. We shall give a general protocol to find this optimal value in the case of stochastic VMC forces, where the situation can be more complicated since the FDT is now sensitive to the QMC intrinsic noise. In our initial tests, the time step is set to 1 fs, a large value, which guarantees a quick and effective exploration of the phase space. Moreover $\gamma_{\textrm{BO}}=0$ and $\Delta_{0}=\delta t$, as we discovered that the most efficient simulation is the one which minimizes the damping in the BO sector. $\Delta_{0}$ is taken as the minimal value which provides a positive definite $\bm{\gamma}^{\textrm{BO}}$. In the top panel of Figure 3, we first observe that in the VMC case, the virial and the primitive kinetic energy estimators deteriorate when the applied input friction is too large ($>0.02\times 10^{-3}$ a.u.). In this situation, the additional QMC noise makes the coupling between the system and the thermostat too large to be fully controlled for such time step values ($\delta t=1$ fs). On the other extreme, at very small $\gamma_{0}$, we see that the presence of a QMC-correlated noise tends to flatten the sharp and deep minimum of the potential autocorrelation time obtained with deterministic CCSD(T) forces (bottom panel). It indicates that, contrary to the deterministic case where there is a clear advantage to set the input friction $\gamma_{0}$ to its optimal value, there is more freedom to choose this parameter in VMC. Indeed, the autocorrelation time divergence shown in the CCSD(T) case for small values of $\gamma_{0}$ disappears with VMC forces. This is due to the implicit low-value cutoff in the $\bm{\gamma}$ matrix provided by the intrinsic QMC noise, once the QMC-force covariance matrix is converted into an effective friction, according to Eqs.(5) and (46). In practice however, the value of $\gamma_{0}$ cannot be too small either, in order to avoid too cold temperatures shown in the middle panel of Fig. 3. Consequently, we need to take the largest $\gamma_{0}$ before the increase of the potential autocorrelation time $\tau_{V}$ due to the soft-modes overdamping. This will also let us recover an acceptable target temperature (see middle panel). Therefore, $\gamma_{0}=1.46\times 10^{-3}$ a.u seems to be a very good compromise between autocorrelation time $\tau_{V}$, effective temperature, and quality of the phase space sampling revealed by the kinetic energy estimators. All subsequent runs will be performed with that value. It is interesting to note that this is optimal for both deterministic and stochastic forces.

Similarly to the tests performed in Sec. II.3 with analytic CCSD(T) forces, we check here the robustness of our novel PIOUD algorithm with respect to the time step $\delta t$ and to the number of quantum replicas $M$ in the presence of noisy QMC forces. In the upper panel of Figure 4, the difference between the virial and the primitive kinetic energy estimators remains very reasonable with increasing values of the time step $\delta t$, even though the time step error is more important in the stochastic case once compared to the deterministic one. The PIOUD propagator exhibits a smaller difference between these two estimators with QMC forces than the PILE algorithm in the deterministic case. The superior performances of PIOUD will allow us to use large time steps $\delta t$.

Finally, we also check the stability of the PIOUD integration scheme with increasing number of quantum replicas $M$. As we can see in the bottom panel of the Figure 4, the difference between the kinetic and the primitive energy estimators is well controlled up to $M=64$ beads in the stochastic case. On the contrary, the PILE propagator already exhibits signs of instability at this value with deterministic forces (see Figure 1). This is a further proof of the robustness of the PIOUD integrator, which is particularly recommended when one wants to perform a PILD simulation with a large number of beads with any force field, deterministic or not.

IV Results

We apply the methodology described previously in Secs. II and III to perform QMC-based LD and PILD simulations of the Zundel ion in the gas phase at low (50 K) and room (300 K) temperature, with and without the NQE. To benchmark our results, we first carried out CMPC and PIOUD calculations using the CCSD(T) PES provided by Huang et alHuang et al. (2005). In that case, the ionic forces are computed as finite differences of the potential energy with an increment $\delta q=10^{-4}$ Bohr in the geometry. The dynamics is propagated with a time step $\delta t=1$ fs until we obtain fully converged atomic distributions. As previously detailed, the input friction ($\gamma_{\textrm{BO}}$ for CMPC, $\gamma_{0}$ for PIOUD) is set to $1.46\times 10^{-3}$ a.u., in order to both ensure a good thermalization of the system and minimize the potential energy autocorrelation time $\tau_{V}$. Consequently, that value also maximizes the diffusion of the nuclei between two subsequent LD iterations, and guarantees an efficient phase space sampling. We use 128 beads at low temperature whereas 32 beads are enough to fully recover NQE at room temperature, which is in agreement with Ref. Spura et al., 2015.

The QMC-LD (CMPC) and QMC-PILD (PIOUD) simulations are performed using the wave function of Ref. Dagrada et al., 2014 as a starting point for the dynamics. Unless otherwise specified, all calculations are carried out by bead-grouping the VMC parameters with $N_{\textrm{groups}}=1$, namely with the same wave function parameters shared by all the replicas at a given iteration. The QMC statistics is the same as the one used in the stability tests (Sec. III.4). The equations of motion are propagated with a large time step $\delta t=1$ fs to explore the phase space as fast as possible without doing too many evaluations of the wave function. We stopped the production run after 20 ps of dynamics, which is the minimum required to obtain converged resultsSit and Marzari (2005). At the end of each simulation, the average temperature, the virial and primitive kinetic energy estimators, the energy fluctuations, and the evolution of the energy gradients with respect to the electronic parameters are checked to ensure the reliability of the simulation. In order to follow the correct BO energy landscape, we require that the error made on the temperature and the kinetic energy must not exceed $5\%$, whereas the values of the energy gradients with respect to the wave function parameters must be lower than 3-4 times their standard deviation.

In Figures 5 and 6 we present the normalized oxygen-oxygen (g${}_{\text{OO}}$) and oxygen-(excess) proton (g${}_{\text{OH}}$) distributions as a function of the inter-oxygen distance. These radial distribution functions are obtained at low and room temperature with and without quantum nuclei. Both the VMC and benchmark CCSD(T) results are shown.

We first observe a very good overall agreement between the radial correlation functions g${}_{\text{OO}}$ and g${}_{\text{OH}}$ obtained with the reference CCSD(T) calculations and our VMC-based simulations. We notice that at 50 K, quantum effects renormalize the height of the g${}_{\text{OO}}$ peak by a factor of 3, while the g${}_{\text{OH}}$ distribution is renormalized by a factor of 2 with respect to its classical counterpart. The NQE have thus a huge effect in the oxygen-oxygen distribution function because the classical system is almost frozen around its zero temperature equilibrium configuration, whereas the ZPE leads to strong quantum fluctuations even at low temperature.

At 300 K, even though NQE are much less significant for the oxygen-oxygen distribution, they are still very important for the oxygen-proton correlation function, where the shapes of the distributions obtained in the classical and in the quantum case are very different. Indeed, the classical g${}_{\text{OH}}$ distribution is asymmetric around the equilibrium distance, whereas the quantum correlation function is symmetric with much longer tails, indicating the possibility of instantaneous proton hoppings by quantum tunneling. On the contrary, NQE are less dramatic at room temperature for the g${}_{\text{OO}}$ radial distribution, as expected from the greater mass of the oxygen atoms which diminishes the thermal de Broglie wavelength.

In the classical case, the intrinsic noise present in the QMC forces tends to spoil the sampling of the instantaneous oxygen-oxygen and oxygen-proton distances, which has to be accurate because the resulting distributions are extremely sharp. In the quantum case, instead, the QMC noise is helpful in improving the quantum delocalization of the nuclei in the desired regions of the phase space. The phase space sampling efficiency seems to be enhanced in the quantum case with respect to the classical one. The quantum results are unexpectedly easier to converge, and they yield radial distributions with reduced error bars compared to their classical counterparts.

Our benchmark system is ideal also to check out the quality of the bead-grouping technique described in Sec. III.2 for the electronic parameters in ab initio VMC-PIOUD simulations. At room temperature the bead-grouping with $N_{\textrm{groups}}=1$ works very well, giving results on the top of the CCSD(T) reference. At low temperature, the agreement between the CCSD(T) reference and the VMC-PIOUD results with $N_{\textrm{groups}}=1$ is still good, although some minor discrepancies appear in the tails of the oxygen-oxygen and oxygen-proton distributions. The strongest bias, though still quantitatively acceptable, is present in the g${}_{\text{OH}}$ function, as this pair distribution is the most sensitive to quantum delocalization effects. By increasing $N_{\textrm{groups}}$ to 16, we improve the peak positions of both g${}_{\text{OH}}$ and g${}_{\text{OO}}$, and the error made in the their tails becomes almost negligible.

This can be simply interpreted by considering the quantum-to-classical isomorphism: hydrogen atoms have a light mass, so the corresponding ring polymers are much more spread than the ones mimicking the quantum nuclei of oxygen. Consequently, the bead-grouping approximation on the optimal electronic parameters $\{\bm{\lambda}_{\textrm{av}}^{(l)}\}_{l=1,\dots,N_{\mathrm{groups}}}$ is more severe for hydrogen than for oxygen. The resulting potential energy landscape is thus affected, and displays a larger curvature around its minimum, due to the energy penalty given by the non fully-optimized wave functions, being the worst for those beads which are the farthest from the centroid. This effect is apparent in the slightly shorter tails of the radial correlation function, since the corresponding ionic configurations are less visited as they have higher energies. Therefore, a compromise must be found by minimizing the total amount of CPU time spent for a simulation and the desired target accuracy on the structural and static properties of the system. While at 50 K an $N_{\textrm{groups}}>1$ should be chosen, at 300 K $N_{\textrm{groups}}=1$ gives very accurate results. Thus, at room temperature we are able to carry out a reliable and very accurate fully-quantum dynamics of the Zundel cation in almost the same CPU time as for classical nuclei.

Finally, it is also interesting to quantify NQE versus thermal effects, by observing the evolution of the radial correlation functions at increasing temperature or when the quantum delocalization of the nuclei is taken into account. For classical particles, there is a clear broadening of the g${}_{\text{OO}}$ and g${}_{\text{OH}}$ distribution with increasing temperature. Indeed, the entropy is increased and the instantaneous oxygen-oxygen and oxygen-proton distances are subject to enhanced fluctuations. This is confirmed by examining Figure 7, where we plot the oxygen-proton distance for the two oxygen sites as a function of the inter-oxygen distance. Dagrada and coworkersDagrada et al. (2014) have shown that QMC techniques are in excellent agreement with the reference CCSD(T) geometries at zero temperature. The minimum energy configuration of the Zundel cation is symmetric (C2 symmetry) with the proton at the middle at the oxygen-oxygen distance. At $T=50$ K, the classical proton remains extremely localized around the C2-symmetry geometric minimum which leads to a very sharp distribution, reported in the upper left panels of Figures 5 and 6. At $T=300$ K, the system has more thermal energy to visit asymmetric configurations with longer oxygen-oxygen distance and the excess proton sitting closer to one oxygen than to the other. This is represented by the two branches in the lower left panel of Figure 7. Even in this situation the symmetric configurations are more often explored than the phase regions forming the two wings.

On the contrary, thermal effects are much less important when NQE are taken into account, as apparent from the right panels of Figures 5, 6 and 7, which show very similar correlation functions for both temperatures. Indeed, quantum fluctuations make the proton able to easily jump to a neighbor oxygen site by quantum tunneling, recovering a more symmetric behavior. This is characterized by longer tails in the g${}_{\text{OH}}$ distribution function compared with the classical ones at $T=300$ K, and by the absence of wings in Fig. 7. Moreover, the 2D distributions are nearly the same at low and room temperature. These conclusions are in agreement with previous studies on this systemSpura et al. (2015); Suzuki et al. (2013). Some calculations performed with deterministic forces obtained from the reference CCSD(T) PES show that thermal effects start to be significant for temperatures greater than $T=900$ K. This confirms that NQE are essential to fully understand the microscopic mechanisms involving the proton dynamics in liquid water at ambient conditions.

Last but not least, we discuss the computational cost of our VMC-(PI)LD simulations. The CPU time required on HPC Marconi at CINECA (2.3 GHz 2 x 18-cores Intel Xeon E5-2697 v4 processors) to carry out a 20 ps trajectory of the Zundel ion using our VMC-(PI)LD approach is about 270k core hours. There is no significant difference between classical and quantum simulations, as long as the path integral dynamics is performed with the bead-grouping approximation with $N_{\textrm{groups}}=1$. The target statistics has been detailed above and is enough to have stable and converged simulations in both electronic and ionic sectors. At the given target statistics, relaxing the bead-grouping approximation in the quantum case implies a total computational time increased by a factor $1\leq N_{\textrm{groups}}\leq M$ with respect to the fastest case of $N_{\textrm{groups}}=1$. In the perspective of applying our novel methodology to larger protonated water clusters, it is worth mentioning that for single-point calculations of equilibrium geometry at zero temperature, it has been estimatedDagrada et al. (2014) that the six-molecule complex represents a crossing point in the relative efficiency between VMC and CCSD(T) methods. By considering the excellent performances of QMC for parallel computations, this technique is already competitive for small systems as far as the elapsed computational time is concerned. Moreover, it is obvious that a further increase in the cluster size would make the QMC approach considerably favored even in terms of total computational demand, thanks to its milder scaling with the system size.

V Conclusions

In this paper, we have extended the previous QMC-driven classical LD formalism (CMPC) to quantum nuclei. We have demonstrated the possibility to deal with the intrinsic QMC noise in the ionic forces evaluation by incorporating this possible source of bias into a PILD formalism. By fulfilling the FDT in a generalized framework, we correct the quantum dynamics for the QMC noise. We are thus able to perform simulations of QMC correlated electrons coupled with quantum nuclei.

From the theoretical side, we work with joint momentum-position coordinates and exploit a Trotter breakup between the physical vibrations, set by the BO surface, and the harmonic ones, introduced by the quantum-to-classical isomorphism. This allows us to separate the dynamics driven by the BO forces from the harmonic one of the free ring polymers coupled with the Langevin thermostat. In our approach, the latter dynamics is propagated exactly, as the resulting stochastic differential equations are linear in both momenta and positions. This step corresponds to an Ornstein-Uhlenbeck process, while in the case of QMC forces, the noise correction is performed in the former BO step. The resulting algorithm is the Path Integral Ornstein Uhlenbeck Dynamics (PIOUD).

As a first application, we took the Zundel ion, a benchmark system, and compared its quantum nuclear dynamics driven by VMC forces with the one evolved according to a deterministic PES, obtained at the accurate CCSD(T) level of theory. In the case of deterministic CCSD(T) forces, we have shown that our novel PIOUD algorithm has properties similar to the PILE scheme as far as the potential autocorrelation time $\tau_{V}$ and the algorithmic diffusion constant $D$ are concerned. Moreover, PIOUD features a larger stability range as a function of both time step $\delta t$ and number of quantum replicas $M$. This is encouraging, since with its improved stability, PIOUD can visit the phase space in a more efficient way, no matter whether the forces are deterministic or not.

We further verified the stability of our approach in the stochastic case, where the PES is computed by VMC during a short nuclear dynamics. Then, we carried out extensive VMC simulations at low (50 K) and room (300 K) temperatures of the Zundel complex with both classical and quantum nuclei, by means of the CMPC and PIOUD algorithms, respectively. The outcome of the calculations in this model system suggests that even at room temperature the NQE are crucial to reach a quantitative description of proton transfer in water and aqueous systems. Indeed, NQE take over thermal effects not only at 50 K but also at 300 K, and contribute to a greater proton mobility by quantum tunneling. A very good agreement with the CCSD(T) reference calculations is found, even in the coarsest (and most effective) bead-grouping approximation, where the same electron VMC parameters are shared by all quantum replicas.

This work proves the potential of QMC methods combined with the PILD approach to perform fully quantum calculations of water-like and aqueous systems. In the perspective of performing such simulations on larger systems, our novel scheme overcomes two major scalability problems. On the one hand, it benefits from the QMC reasonable scaling with respect to the number of electrons contrary to other advanced electronic structure methods. On the other hand, using the bead-grouping technique, we are able to carry out fully ab initio dynamics of the Zundel ion with almost the same computational cost as for classical nuclei, without deteriorating NQE. Within this approximation, the quantum simulations are paradoxically more efficient than their classical counterparts. Indeed, the nuclear observables statistics is improved, because the phase space is more efficiently explored, thanks to the quantum fluctuations included in the framework, while the QMC statistical fluctuations of the electronic part are not detrimental, as they are averaged over the beads. Consequently, this work paves the way to study, at the theoretical level, proton transfer mechanisms in more complex water clusters, starting, for instance, with the protonated water hexamer $\text{H}_{13}\text{O}_{6}^{+}$, and the larger $\text{H}_{43}\text{O}_{21}^{+}$ cluster, where the excess proton localization is still under debate. Similar ideas developed in this work for the application of the molecular dynamics to describe NQE, can be used within the recently proposed ”accelerated molecular dynamics” based on first order Langevin equation.Mazzola and Sorella (2017)

VI Associated content

VI.1 Supporting Information

Comparison between the PIOUD and PILE algorithms; average simulation temperature, algorithmic diffusion constant, and the potential autocorrelation time as a function of both number of beads and time step; two figures illustrating their behavior for the simulation of the Zundel ion at 300 K. This material is available free of charge via the Internet at https://pubs.acs.org/doi/suppl/10.1021/acs.jctc.7b00017.

Acknowledgements

We thank Joel M. Bowman and Xinchuan Huang for providing us with the coupled cluster potential energy surface of $\text{H}_{5}\text{O}_{2}^{+}$. F.M. thanks Mario Dagrada for sharing with us the optimal $\text{H}_{5}\text{O}_{2}^{+}$ wave function at zero temperature, and for helpful discussions. We also acknowledge computational resources provided by the PRACE preparatory access number 2010PA2475 and PRACE projects 2015133179 and 2016133322. This work was granted access to the HPC resources of TGCC and IDRIS under the allocation 2016-096493 made by GENCI.

Appendix A Noise correlators in the CMPC algorithm

We report the variance-covariance properties which define the time integrated noises $\bm{\tilde{\eta}}$ and $\bm{\tilde{\tilde{\eta}}}$ for momenta and positions, respectively. Their variance-covariance can be computed by imposing the FDT (5) to be fulfilled, under the hypothesis that the $\bm{\alpha}$ matrix is $\mathbf{q}$-time independent (as we have done in Sec. II.1), and by exploiting that $\left[\bm{\alpha},\bm{\gamma}\right]=0$. One finds:

[TABLE]

while the mean of $\bm{\tilde{\eta}}$ and $\bm{{\tilde{\tilde{\eta}}}}$ is zero.

Appendix B QMC noise correction in the CMPC algorithm

We provide the noise correction which needs to be applied if one wants to include the QMC noise into the Langevin dynamics without inducing a bias on the final target temperature. Analogously to the equation (48), the integrated noise can be written

[TABLE]

Therefore, still assuming that the external noise is independent from the QMC noise, we obtain very similar relations to Eq.(50):

[TABLE]

where the explicit expressions of the $\bm{\tilde{\alpha}}$ matrix components have been given in the main text in Eq. (A). By simply computing the square root of this $2\times 2$ matrix which is positive definite by construction, we are able to evaluate explicitly the exact integrated external noise to propagate the dynamics.

Appendix C Quantum integration scheme

We detail in this section the mathematical derivation of our novel integrator described in Section II.2 and the explicit formulas of the $\bm{\Lambda}$, $\bm{\Gamma}$ and $\bm{\Theta}$ matrices. For the following algebra, it is useful to evaluate the inverse of the matrix $\bm{\hat{\gamma}}$

[TABLE]

where we replaced $\bm{\gamma}^{\textrm{harm}}$ of Eq. 25 by $\bm{\gamma}$, for the sake of readability. In order to solve the differential system in Eq. (10) for a generic time step $\delta t$ in the quantum case, we need to exponentiate the matrix $\bm{\hat{\gamma}}$. Using the fundamental assumption $\left[\mathbf{K},\bm{\gamma}\right]=0$ previously justified in the paper, we consider each common eigenvector of $\mathbf{K}$ and $\bm{\gamma}$ that can correspond to a joint momentum and coordinate mode. In this two-fold basis, the matrix $\bm{\hat{\gamma}}$ is a simple $2\times 2$ block matrix where as $K$ and $\gamma$ are just numbers. The block matrix $\bm{\hat{\gamma}}$ can be more conveniently rewritten in terms of the Pauli matrices $\bm{\sigma}_{x}$,$\bm{\sigma}_{y}$,$\bm{\sigma}_{z}$:

[TABLE]

Then, the exponentiation can be straightforwardly obtained, by using standard Pauli matrices algebra

[TABLE]

where $|x|=\sqrt{|\gamma^{2}/4-K|}$. Recombining Eq.(54) and Eq.(56) one obtains:

[TABLE]

for $\gamma^{2}>4K$, while for $\gamma^{2}\leq 4K$:

[TABLE]

Finally, for the case $\gamma^{2}>4K$, we obtain:

[TABLE]

while for $\gamma^{2}\leq 4K$, we have:

[TABLE]

To evaluate the time-integrated noise covariance matrix $\bm{\tilde{\alpha}}$ in the quantum case, defined by Eq. (31), namely $\bm{\tilde{\alpha}}\equiv\langle\mathbf{E}_{\textrm{int}}^{T}\mathbf{E}_{\textrm{int}}\rangle$, it is useful to introduce the following integral:

[TABLE]

By using Eq. (56), it is straightforward to carry out the integrals to obtain

[TABLE]

for $\gamma^{2}\geq 4K$, whereas for $\gamma^{2}\leq 4K$ we obtain a very similar formula where $|x|$ is replaced by $i|x|$:

[TABLE]

$\bm{\hat{\alpha}}$ has been defined in Eq. (32). In order to conclude the analytic derivation of the quantum integration scheme, we also give the explicit expression of the two matrix products appearing in the Eq. (62) and Eq. (63):

[TABLE]