Path-integral dynamics of water using curvilinear centroids

TL;DR

This paper introduces a path-integral dynamics method for water using curvilinear centroid coordinates, improving quantum spectral simulations and avoiding common CMD issues, with potential applications to other systems.

Contribution

It develops a novel curvilinear centroid molecular dynamics approach that enhances quantum simulations of water and potentially other molecular systems.

Findings

QCMD spectra closely match exact quantum spectra across conditions

QCMD avoids the CMD curvature problem in spectral predictions

QCMD results are consistent with other quantum simulation methods

Abstract

We develop a path-integral dynamics method for water that resembles centroid molecular dynamics (CMD), except that the centroids are averages of curvilinear, rather than cartesian, bead coordinates. The curvilinear coordinates are used explicitly only when computing the potential of mean force, the components of which are re-expressed in terms of cartesian 'quasi-centroids' (so-called because they are close to the cartesian centroids). Cartesian equations of motion are obtained by making small approximations to the quantum Boltzmann distribution. Simulations of the infrared spectra of various water models over 150-600 K show these approximations to be justified: for a two-dimensional OH-bond model, the quasi-centroid molecular dynamics (QCMD) spectra lie close to the exact quantum spectra, and almost on top of the Matsubara dynamics spectra; for gas-phase water, the QCMD spectra are close to the exact quantum spectra; for liquid water and ice (using the q-TIP4P/F surface), the QCMD spectra are close to the CMD spectra at 600 K, and line up with the results of thermostatted ring-polymer molecular dynamics and approximate quantum calculations at 300 and 150 K. The QCMD spectra show no sign of the CMD 'curvature problem' (of erroneous red shifts and broadening). In the liquid and ice simulations, the potential of mean force was evaluated on the fly by generalising an adiabatic CMD algorithm to curvilinear coordinates; the full limit of adiabatic separation needed to be taken, which made the QCMD calculations 8 times more expensive than partially adiabatic CMD at 300 K, and 32 times at 150 K (and the intensities may still not be converged at this temperature). The QCMD method is probably generalisable to many other systems, provided collective bead-coordinates can be identified that yield compact mean-field ring-polymer distributions.

Equations

$$\mathbf{Q}=\frac{1}{N}\sum _{i=1}^N{\mathbf{q}}_i$$

$$\dot{\mathbf{Q}}=$$

$$\dot{\mathbf{P}}=$$

$$F(\mathbf{Q})=$$

$$Z_0(\mathbf{Q})=$$

$$W(\mathbf{q})$$

$$U(\mathbf{q})$$

$$S(\mathbf{q})$$

$$-\frac{\partial F(\mathbf{Q})}{\partial \mathbf{Q}}=-{⟨\frac{\partial U(\mathbf{q})}{\partial \mathbf{Q}}⟩}_{\mathbf{Q}}$$

$${⟨\dots ⟩}_{\mathbf{Q}}$$

$$C(t)=$$

$$H(\mathbf{P},\mathbf{Q})=\frac{{\mathbf{P}}^2}{2m}+F(\mathbf{Q})$$

$$R$$

$$r_i$$

$${\overline{Q}}_x$$

$${\overline{Q}}_y$$

$$\dot{\overline{\mathbf{Q}}}$$

$$\dot{\overline{\mathbf{P}}}$$

$$\overline{F}(R)$$

$${\overline{Z}}_0(R)$$

$$-\frac{\partial \overline{F}(R)}{\partial \overline{\mathbf{Q}}}\simeq -{⟨\frac{\partial U(\mathbf{q})}{\partial \overline{\mathbf{Q}}}⟩}_R$$

$${⟨\dots ⟩}_R$$

$$-\frac{\partial U(\mathbf{q})}{\partial \overline{\mathbf{Q}}}$$

$$f_R(\mathbf{q})$$

$$\overline{C}(t)$$

$$\overline{Z}$$

$$\overline{H}(\overline{\mathbf{P}},\overline{\mathbf{Q}})={\overline{\mathbf{P}}}^2/2m+f(\overline{\mathbf{Q}})$$

$$R_1$$

$$R_2$$

$$\Theta$$

$${\dot{\overline{\mathbf{Q}}}}^{(\alpha )}$$

$${\dot{\overline{\mathbf{P}}}}^{(\alpha )}$$

$$-\frac{\partial \overline{F}(\mathbf{\xi })}{\partial {\overline{\mathbf{Q}}}^{(\alpha )}}\simeq -{⟨\frac{\partial U(\mathbf{q})}{\partial {\overline{\mathbf{Q}}}^{(\alpha )}}⟩}_{\mathbf{\xi }}.$$

$${⟨\dots ⟩}_{\mathbf{\xi }}$$

Full text

I Introduction

A major challenge in molecular simulation is to include nuclear quantum effects in the dynamics of liquid water, ice, and related systems. One way to do this is to make local approximations around one water molecule, then to solve the Schrödinger equation in a reduced space.skinner2009 ; skinner2013 ; lmon1 ; lmon2 Another way, arguably more general, because it does not require a local approximation, is to use path-integrals to describe the quantum Boltzmann statistics exactly, and classical mechanics to approximate the dynamics. billjpc ; liubill ; liurev ; shig ; jens1 ; xantheas ; liuwat ; bonella ; vollume ; van1 ; van2 ; jerin1 ; jerin2 ; caov2 ; craig1 ; craig2 ; tfmiller1 ; annurev ; trpmd1 ; jorsca ; tim ; tom1 ; tom2 ; paes1 ; paes2 This type of approach assumes that real-time coherences mainly average out in condensed-phase systems, leaving most of the quantum effects in the statistics.

The oldest such ‘quantum statistics + classical dynamics’ approach is based on the classical Wigner, or LSC-IVR, approximation.billjpc ; liubill ; liurev ; shig ; jens1 This approach is powerful, and has been used to calculate infraredxantheas and recently Raman spectraliuwat for liquid water. However, LSC-IVR has the drawback that the classical trajectories do not conserve the quantum Boltzmann distribution, which has the effect of damping spectral lineshapes.

More recently, approaches such as centroid molecular dynamics (CMD)caov2 and (thermostatted) ring-polymer molecular dynamics ((T)RPMD)craig1 ; craig2 ; tfmiller1 ; annurev ; trpmd1 ; jorsca ; tim ; tom1 ; tom2 have been developed, which redress this problem by using trajectories that do conserve the quantum Boltzmann distribution. These approaches are very practical, since they can be run as extended classical MD simulations. But they also have drawbacks of their own.

The CMD method has proved to be useful for simulating water dynamics, most notably in recent calculations on the MB-pol ab initio water potential energy surface.paes1 ; paes2 The method works by generating classical trajectories on the free-energy surface obtained by mean-field averaging the quantum Boltzmann distribution around the centroids (centres-of-mass) of the ring-polymers (imaginary-time Feynman pathschanwol ; pariram ; ceperley ). However, CMD breaks down at low temperatures, when the mean-field ring-polymer distributions delocalise.trpmd2 ; marx09 ; marx10 When applied to gas-phase water, the CMD ring-polymer distribution spreads in a crescent around the curve of the rotating OH-bond, so that the centroid lies approximately at its focal point, which flattens the potential, and leads to an artificially red-shifted and broadened spectral lineshape. This ‘curvature problem’ is less severe in the condensed than the gas phase,marx10 since there is less spreading around librations than rotations, but it is still serious enough to give, e.g., a 100 cm*-1* red shift and considerable broadening to the OH-stretch band of ice at 150 K.trpmd2

The RPMD method does not suffer from the CMD curvature problem, since it involves no mean-field averaging, propagating instead all the ring-polymer ‘beads’ (imaginary-time replicas) as separate particles. However, the dynamics of the modes orthogonal to the centroid are fictitious, because of the ring-polymer spring-forces, and hence need to be damped by a thermostat to avoid corrupting the spectrum.trpmd1 ; trpmd2 ; trpmd3 This means that thermostatted (T)RPMD is an extremely efficient method for generating spectral line positions, but that it predicts lineshapes that are often too broad.

It would therefore be useful to develop a path-integral method which generates dynamics that is quantum Boltzmann conserving, but which manages to avoid both the curvature problem of CMD, and the fictitious springs of RPMD. In refs. marx09, ; marx10, , Marx and co-workers suggest that such a method could be developed by modifying CMD, such that the centroids are taken as the averages of curvilinear rather than cartesian bead coordinates. This curvilinear centroid would lie at the centre of the ring-polymer distribution, rather than at its approximate focal point, which would stop the free energy from artificially flattening. To our knowledge, this suggestion has not yet been pursued.

More recently, the case for curvilinear centroids has been strengthened further, by the development of ‘Matsubara dynamics’, which is the quantum-Boltzmann-conserving classical dynamics that results when the exact quantum dynamics is mean-field averaged to make it smooth and continuous in imaginary time.mats ; mf-mats Matsubara dynamics is impossible to compute exactly (except in toy systems), owing to a phase problem, but it can be computed approximately.jens2 ; planets ; mats-rpmd In fact CMD and TRPMD (both first obtained heuristically) are mean-field and short-time approximations to Matsubara dynamics.mats-rpmd By considering a two-dimensional model, it was shown that the CMD curvature problem in the OH stretch can be cured if just a few modes are released from the centroid-averaged mean-field.mf-mats The resulting Matsubara stretch bands are very close to the exact quantum results, the only significant error being a small, temperature-independent blue shift, which also affects CMD (when the temperature is high enough to avoid the curvature problem) and TRPMD. It was also found that the spreading of the CMD ring-polymer distribution around the rotational angle is exacerbated by artificial instantons, which form at small bond-lengths because of the cartesian centroid constraint.mf-mats Curvilinear centroids would avoid such artefacts, and should therefore give more compact distributions around the rotational angle, and thus a better approximation to the Matsubara-dynamics spectrum.

In this article, we investigate these likely advantages, by developing a method, called ‘quasi-centroid molecular dynamics’ (QCMD), which applies curvilinear centroid dynamics to water. We have not tried formally to derive QCMD from Matsubara dynamics—such a derivation is almost certainly possible, but we thought it best first to test QCMD numerically.

Before going further, however, we mention some of the concerns that have probably put people off using curvilinear centroids to date. First, curvilinear hamiltonians are notoriously difficult to derive and work with,miller-ham ; tennyson-ham ; truhlar-ham especially when used for path-integrals (which exhibit strange pathologies when formulated in curvilinear coordinateskleinert ; marx-curv ). Second, in using non-cartesian centroids, one runs the risk that components of the ring-polymer springs could survive the mean-field averaging, and thus corrupt the forces. Third, except in the high-temperature limit, the overall curvilinear centroid (which we will refer to as the ‘quasi-centroid’) does not coincide with the overall cartesian centroid (which lies at the bead centre-of-mass); as a result, the quasi-centroid does not give the exact estimator for a linear operator. Our response to these concerns is to formulate the mean-field dynamics in as cartesian a way as possible, by making approximations to the quantum Boltzmann distribution which are expected to be minor, provided it is sufficiently compact that the quasi-centroids are close to the centroids. The numerical results in Sec. II, III, and V will show that these assumptions appear to be justified.

The article is structured as follows. Section II develops the QCMD method for a two-dimensional OH-stretch model, for which the curvilinear centroid is simply the average of the bead radii. The resulting QCMD spectra are found to lie close to the exact quantum spectra, and almost on top of the Matsubara spectra. Section III extends these results to gas-phase water, by using bond-angle centroids, obtaining similarly close agreement with the exact quantum results. Section IV shows how to compute the mean-field quasi-centroid forces on-the-fly, using a modification of the adiabatic CMD (ACMD) algorithm,caov4 ; hone1 ; hone2 with SHAKE and RATTLEryckaerts ; andersen ; tuckerman used to impose curvilinear constraints. Section V reports applications of QCMD to liquid water and ice, using the bond-angle centroids of Sec. III to describe the internal motion of the monomers, and an Eckart-like rotational frameeckart ; bright-wilson to describe their translational and librational motion. Section VI concludes the article with suggestions of how to generalise QCMD to systems other than water.

II Curvilinear centroids in a two-dimensional model

To develop the key ideas behind QCMD, we first consider a familiar two-dimensional ‘champagne-bottle’ model (details given in Appendix A) of a vibrating-rotating OH bond, with a linear dipole moment operator.

II.1 Standard CMD treatment

In CMD, one considers the motion of the centroid, which is the centre of mass of the ring-polymer.caov2 For a particle of mass $m$ moving in two dimensions, represented by $N$ ring-polymer beads, the centroid position ${\bf Q}\equiv(Q_{x},Q_{y})$ is

[TABLE]

where ${\bf q}_{i}\equiv(q_{ix},q_{iy})$ are the cartesian coordinates of the $i$-th ring-polymer bead. The CMD equations of motion are

[TABLE]

where ${\bf P}$ is the centroid momentum, and

[TABLE]

is the free energy obtained fromreduced

[TABLE]

in which

[TABLE]

Since the spring potential $S({\bf q})$ is independent of ${\bf Q}$, it follows that

[TABLE]

where

[TABLE]

which is useful in practical calculations, since the right-hand side can be evaluated using methods such as adiabatic centroid molecular dynamics (ACMD),caov4 ; hone1 ; hone2 as discussed in Sec. IV.

The CMD approximation to the infrared spectrum is obtained from the Fourier transform of

[TABLE]

with

[TABLE]

In the model of Appendix A, we use a linear dipole moment $\bm{\mu}=\mathbf{Q}$, and multiply $C(t)$ by a time-window function before taking the Fourier transform.

Figure 1 illustrates the performance of CMD for the two-dimensional OH-bond model. As is well known, CMD works well for this model at relatively high temperatures, but the ‘curvature problem’ artificially red-shifts and broadens the OH stretch at lower temperatures, where the centroid-constrained ring polymers spread around the potential curve.marx09 ; marx10 In ref. 38, it was shown that this behaviour results from artificial instantons minimising the ring-polymer potential energy at small bond-lengths. The curvature problem worsens noticeably once the temperature is low enough for the CMD trajectories to encounter these instantons, as happens at about $200~{}$\mathrm{K}$$ in the OH model (see Fig. 1 showing the spectra and Fig. 3 showing an instanton-forming CMD trajectory). That the instantons minimise the ring-polymer energy is an artefact of applying a cartesian constraint to a curved potential: by stretching the polymers round the curve, one can push them outwards (away from the repulsive wall) without affecting the constraint.

II.2 Quasi-centroid molecular dynamics

We therefore propose constraining the centroid

[TABLE]

of the ring-polymer radial coordinates

[TABLE]

This constraint makes it impossible for the polymers to lower their energy by stretching and moving outwards into instantons, and ensures that $R$ describes the centre of the ring-polymer distribution rather than its approximate focal point.instantons We also define an analogous polar angle $\Theta$, which runs from [math] to $2\pi$ such that the cartesian coordinates

[TABLE]

cover the whole two-dimensional space. Note that the dependence of $\Theta$ on the bead coordinates need not be specified, because the circular symmetry of the model means that the potential of mean force is independent of $\Theta$.

The non-linear relation between the cartesian and polar bead coordinates means that $\overline{\bf Q}\neq{\bf Q}$. However, it is easy to show (by expanding $R$ in terms of ring-polymer normal modes) that $\overline{\bf Q}\simeq{\bf Q}$ if the distribution is reasonably compact, and that $\overline{\bf Q}\to{\bf Q}$ in the high temperature limit. We will therefore refer to $\overline{\bf Q}$ as the position of the ‘quasi-centroid’.

In quasi-centroid molecular dynamics (QCMD), we use cartesian equations of motion which resemble those of CMD, except that the mean-field averages are taken around the quasi-centroid rather than the centroid:

[TABLE]

where

[TABLE]

with

[TABLE]

In writing out these equations we have assumed that it is a good approximation to mean-field average the exact quantum dynamics about the quasi-centroid phase space $(\overline{\bf P},\overline{\bf Q})$, and that this mean-field-averaged dynamics is classical. The first assumption is expected to hold if the distribution is compact, and the second one if the mean-field-averaging gives a phaseless Matsubara dynamics (which we have not derived, but which seems likely). In writing Eq. (15b), we have also assumed that $\overline{\bf P}$ can be approximated to be purely cartesian, such that it contributes a factor of $\exp(-\beta\overline{\bf P}^{2}\!/2m)$ to the quantum Boltzmann distribution. Appendix B shows that this last assumption can be expected to hold if the ring-polymer distribution is sufficiently compact.

We further simplify the calculation of the force by making the approximation

[TABLE]

where

[TABLE]

which is equivalent to assuming that the polymer-spring forces do not survive the mean-field averaging.Jacob Unlike its counterpart in CMD, Eq. (18) is not exact, since $R$ contains components of the ring-polymer normal modes orthogonal to the centroid. However, the size of this spring-force will depend on the extent to which the ring-polymer distribution spreads and contracts as a function of $R$, and is thus expected to be small if the distribution remains compact. The advantage of Eq. (18) is that it will allow us later (Sec. IV) to evaluate the force on-the-fly using a simple algorithm. The cartesian forces of Eq. (18) are obtained from the forces on $R$ using

[TABLE]

with

[TABLE]

Having propagated the quasi-centroid trajectory $\overline{Q}(t)$, the spectrum is then obtained from the time-correlation function

[TABLE]

with

[TABLE]

and

[TABLE]

where $f(\overline{\mathbf{Q}})\simeq\overline{F}(\overline{\mathbf{Q}})$ is defined to be the free energy obtained by doing work with the approximate force of Eq. (18).

In a QCMD simulation, one propagates trajectories using Eqs. (15) and (18), applying a standard thermostat in order to sample the distribution $\exp(-\beta\overline{H}(\overline{\bf P},\overline{\bf Q}))$. As already mentioned, this distribution involves approximations to the exact quantum Boltzmann distribution which would not have been made in the analogous CMD calculation. In addition, CMD has the advantage of using the exact estimator for the linear dipole moment operator ${\bm{\mu}}={\bf Q}$, whereas QCMD uses the approximation ${\bm{\mu}}\simeq\overline{\bf Q}$. We thus expect QCMD to be less accurate than CMD when the temperature is high enough to neglect the curvature problem, but more accurate at lower temperatures, where the better treatment by QCMD of the dynamics will far outweigh the errors in the Boltzmann distribution which are expected to be small (provided the ring-polymer distribution is compact—see Appendix B). To be safe, we can monitor the errors in the distribution numerically, by comparing QCMD static properties with results from standard path-integral molecular dynamics (PIMD) calculations.

We tested QCMD on the two-dimensional OH-bond model of Appendix A. The forces in Eq. (21) were evaluated at 128 grid-points along the $R$ coordinate, using standard PIMDtuckerman with the quasi-centroid constraint imposed using SHAKE and RATTLE.ryckaerts ; andersen ; tuckerman The resulting spectra (Fig. 1) are in excellent agreement with the exact quantum and Matsubara results, at all temperatures tested, with no sign of a curvature problem. Figure 2 makes a closer comparison of the QCMD OH-stretch bands with the mean-field Matsubara results of ref. 38. We see that the QCMD OH-stretch peaks have small ($\sim\!10~{}\mathrm{cm}^{-1}$) blue shifts with respect to the Matsubara peaks. The shifts change slightly with temperature, ranging from $14~{}\mathrm{cm}^{-1}$ at $800~{}\mathrm{K}$ to $6~{}\mathrm{cm}^{-1}$ at $200~{}\mathrm{K}$. This variation appears to be the result of small differences in lineshape, rather than an underlying curvature problem; part of these differences may be the result of convergence artefacts in the Matsubara results.

As expected, the QCMD distributions contain no artificial instantons at the inner turning points (see Fig. 3). Detailed comparison with CMD trajectories is difficult (since the dynamics are different), but it appears that the QCMD distributions are broader than the CMD distributions at the outer turning points. This property may perhaps be responsible for the small blue shifts in comparison with the Matsubara results (see Fig. 2)—although if this were case, one would expect a systematic temperature dependence. Overall, the widths of the QCMD distributions vary little over a vibrational period, confirming that dynamics of the quasi-centroid is coupled only weakly to the internal degrees of freedom, so that the neglect of the ring-polymer springs in Eq. (18) is justified.

Given the excellent agreement of the QCMD and Matsubara spectra, we expect that the various QCMD static approximations mentioned above are small. This is confirmed by comparing QCMD static averages with the results of exact quantum calculations. We found that the QCMD expectation values of the bond-length $R$ were identical to within sampling error to the exact values, at all temperatures tested, and that the largest error in the $t=0$ limit of the dipole autocorrelation function was 2$\%$ at 200 K. Note that this error does not affect the calculated spectra, since these were obtained from the dipole-derivative autocorrelation function, for which QCMD gives the exact $t=0$ limit for a linear dipole.

We attempted to ‘break’ the QCMD calculations by pushing the temperature down to 50 K. However, even at this low temperature, the ring-polymer distribution around the centroid remained compact, and the agreement of QCMD with the exact quantum spectrum was as close as that shown at higher temperatures in Fig. 1. This suggests that QCMD continues to gives a good approximation to Matsubara dynamics down to 50 K (where we were unable to compute the Matsubara spectrum because the phase problem was too severe).

In summary, QCMD behaves extremely well for the two-dimensional OH model. It eliminates the curvature problem, giving spectra that lie almost on top of the Matsubara results at all temperatures tested. The errors from the ‘extra’ approximations (that $\overline{\bf P}$ is cartesian and that the polymer springs can be neglected in the force) are found to be insignificant.

III Extension to gas-phase water

We now test whether these advantages continue when QCMD is extended to treat gas-phase water. By analogy with Sec. II, we choose a set of curvilinear centroid coordinates that can be expected to give compact mean-field distributions, namely the bond-angle centroids

[TABLE]

in which $\mathbf{q}^{\scriptscriptstyle{(\mathrm{OH}_{1,2})}}_{i}={\bf q}^{\scriptscriptstyle{(\mathrm{H}_{1,2})}}_{i}\!-{\bf q}^{\scriptscriptstyle{(\mathrm{O})}}_{i}$​ are the displacement vectors between the oxygen and hydrogen atoms in the $i$-th water molecule replica ($i=1\ldots N$). A further six coordinates would be required to orient the mean-field-averaged water molecule and locate its centre-of-mass, but, like $\Theta$ in Sec. II, these external coordinates need not be specified, since no force or torque acts on them.

As in Sec. II, we carry out dynamics using equations of motion which have the same form as CMD, namely

[TABLE]

where ${\overline{\bf Q}}^{(\alpha)}$ locate the quasi-centroids of atoms $\alpha=({\rm H}_{1},{\rm H}_{2},{\rm O})$ in the mean-field averaged water molecule. As in Sec. II, the curvilinear coordinates enter only through the force, which, by analogy with Eq. (18), is approximated as

[TABLE]

where

[TABLE]

with $\bm{\xi}\equiv\left(\xi_{1},\xi_{2},\xi_{3}\right)\equiv\left(R_{1},R_{2},\Theta\right)$, and ${\bf q}\equiv\left\{{\bf q}_{i}^{(\alpha)}\right\}$.

The right-hand side of Eq. (27) is computed using

[TABLE]

with $(k,j)=(1,2)$ or $(2,1)$, and

[TABLE]

where $\bm{\nabla}_{\!\mathrm{H}_{k}\!}V(\mathbf{q}_{i})$ is the force on the $i$-th bead of the $k$-th hydrogen. The infrared spectrum is obtained from

[TABLE]

with

[TABLE]

where $f(\overline{\mathbf{Q}})\simeq\overline{F}(\overline{\mathbf{Q}})$ is the free energy obtained by doing work with the approximate force of Eq. (27).

We used the above framework to simulate QCMD infrared spectra for gas-phase water using the Partridge-Schwenke potential surfacepspes and dipole moment.psdms As in Sec. II, we damped the time-correlation functions before calculating their Fourier transform, this time using a Hann windownr77 of 750 $\mathrm{f}\mathrm{s}$. The quasi-centroid forces in Eq. (30) were calculated on a $64\times 64\times 64$ grid, using standard PIMDtuckerman with SHAKE/RATTLE to impose the quasi-centroid constraints,ryckaerts ; andersen ; tuckerman and were then fitted to cubic splines.nr77 At 600 and 300 K, the grids ranged from 1.50 to 2.50 a.u. in $R_{1,2}$ and from $85^{\circ}$ to $130^{\circ}$ in $\Theta$; at 150 K this range was reduced to 1.65 – 2.05 a.u. and $90^{\circ}$ – $126^{\circ}$. We also calculated the exact quantum spectrum (subject to the same damping) using the DVR3D package of Tennyson and co-workers,dvr3d and the corresponding CMD and TRPMD spectra using standard PIMD techniques. A comparison with Matsubara dynamics was not possible because the Matsubara sign problem was too severe.

From Fig. 4, we see that QCMD works extremely well for gas-phase water. The overall agreement with the exact quantum results is excellent, even at 150 K, where there is again no sign of a curvature problem. There are small differences between the QCMD and quantum spectra: the QCMD bend ($\sim\!1600~{}\mathrm{cm}^{-1}$) and stretch ($\sim\!3800~{}\mathrm{cm}^{-1}$) have different internal structures from those in the exact spectrum, and the stretch is blue-shifted by about $60~{}\mathrm{cm}^{-1}$. The main cause of these differences is most likely the neglect of real-time coherence by QCMD. Note that the positions of the TRPMD bands are very close to those of the QCMD bands, implying that, as expected, TRPMD gives a good description of the short-time dynamics of gas-phase water.

As in the two-dimensional model, the gas-phase water QCMD distributions remain compact at inner turning points because there are no quasi-centroid-constrained instantons. Figure 5 compares the distributions at inner turning points of an ‘extreme’ QCMD and CMD trajectory at 150 K. As in the two-dimensional case, the QCMD distribution remains localised whereas the CMD distribution spreads around a centroid-constrained instanton. This artefact distorts the approximate axial symmetry of the distribution around the OH bond into a banana shape.

Comparisons of the QCMD average OH bond lengths with those obtained using standard PIMD give tiny differences ($0.1\%$ at 150 K), indicating that, as for the two-dimensional model, the approximate distribution in Eq. (32) is very close to the exact quantum Boltzmann distribution. Note that, since the Partridge-Schwenke dipole moment is non-linear,psdms neither CMD nor QCMD give the exact quantum Kubo time-correlation function in the $t\to 0$ limit,hone2 ; paes2008 and that both methods neglect the Matsubara fluctuation terms necessary to completely describe the dipole-moment dynamics at low temperatures.mats ; mats-rpmd ; mf-mats

IV Implementing QCMD using an adiabatic propagator

In Sec. II and III, we applied QCMD to systems that are small enough for the mean-field force on the quasi-centroid to be pre-calculated on a grid. To apply QCMD to larger systems, we adapt the adiabatic CMD (ACMD) algorithm, used in previous workcaov4 ; hone1 ; hone2 to calculate the mean-field centroid force on-the-fly.

In ACMD, the ring polymer is represented in terms of the normal modes that diagonalise the ring-polymer spring potential. These are Fourier modes, such that for each cartesian degree of freedom the lowest frequency mode $s_{0}$ is (proportional to) the centroid.coalson The mass associated with each non-centroid mode $s_{n}$ is then scaled by a factor of $(\beta\hbar\omega_{n}/\gamma N)^{2}$, where $\omega_{n}$ is the associated normal-mode frequency, and a local thermostat is applied to its dynamics. The limit $\gamma\to\infty$ gives a perfect adiabatic separation between the motion of the centroid and the thermostatted normal modes, and is thus equivalent to CMD.hone2 In practice, $\gamma$ is treated as a convergence parameter, with e.g. $\gamma=64$ giving a good approximation to CMD in gas-phase water at 150 K.

It is straightforward to modify the ACMD algorithm to an AQCMD algorithm that adiabatically separates the motion of the quasi-centroid from that of all the other degrees of freedom. In the AQCMD algorithm, a simulation is treated as two systems evolving in parallel, one of quasi-centroids with phase-space coordinates $(\overline{\mathbf{P}},\overline{\mathbf{Q}})$, and one of ring-polymers with coordinates $(\mathbf{p},\mathbf{q})$. The dynamics is propagated using a modification of the integrator due to Leimkuhler and Matthews with OBABO operator splitting.leimkuhler A propagation step begins with the thermostatting of the momenta:

1. O1.

Propagate the quasi-centroid momenta $\overline{\mathbf{P}}$ for half a time-step ($\Delta t/2$) under the action of a Langevin thermostat. 2. O2.

Propagate the bead momenta ${\bf p}$ for $\Delta t/2$ under the path-integral Langevin equation (PILE) thermostat,pile followed by RATTLEandersen to constrain the quasi-centroid components.

The momenta are then updated according to the forces acting on the system:

1. B1.

Propagate $\overline{\bf P}$ for $\Delta t/2$ under the quasi-centroid forces. 2. B2.

Propagate $\mathbf{p}$ for $\Delta t/2$ under the forces derived from the ring-polymer potential in Eq. (5), followed by RATTLE.

This is followed by a position update:

1. A1.

Propagate the quasi-centroid positions $\overline{\mathbf{Q}}$ for a full time-step $\Delta t$ according to the current values of the momenta $\overline{\mathbf{P}}$. 2. A2.

Propagate the bead positions $\mathbf{q}$ for $\Delta t$ according to the current values of $\mathbf{p}$, followed by SHAKE,ryckaerts which constrains the ring-polymer geometry to be consistent with the quasi-centroid configuration $\overline{\mathbf{Q}}$ obtained at the end of step A1.

The propagation step is concluded by executing B1–B2 using the forces evaluated at the updated positions, followed by O1–O2. The mass-scaling used in the propagation of the ring-polymers is the same as in ACMD, except that the mass of the centroid is also scaled, by $\gamma^{-2}$. A useful feature of this algorithm is that it requires explicit forces only on the quasi-centroid (B1) and on the cartesian bead coordinates (B2), but does not require explicit forces on the degrees of freedom orthogonal to the quasi-centroid. Note that the AQCMD algorithm reduces to the ACMD algorithm in the special case that the quasi-centroid is taken to be the centroid (at which point the constraints become trivial to apply).

We tested the convergence with respect to $\gamma$ by applying the AQCMD algorithm to gas-phase water at 150 K. A value of $\gamma=64$ was sufficient to reproduce the spectrum obtained by interpolating the force on a grid. The adiabatic simulation required a time-step of $\Delta t=0.1/\gamma~{}$\mathrm{f}\mathrm{s}$$, i.e. $\gamma$ times smaller than what would have been used to propagate the dynamics without the mass scaling. We note that the same values of $\gamma$ and $\Delta t$ are required to converge ACMD at this temperature, and hence that the AQCMD algorithm uses approximately the same amount of CPU time as the ACMD algorithm, since most of the time is spent evaluating the forces on the ring-polymer beads.

In condensed phase simulations, converged ACMD results can be obtained at significantly lower adiabatic separation,hone2 ; mano-pacmd provided the strength of the thermostat acting on the non-centroid modes is carefully tuned.trpmd2-sup This approach is known as partially-adiabatic CMD (PACMD). We have not yet attempted a similar partially-adiabatic implementation of QCMD and view it as an important step in future methodological development.

V Liquid water and ice

Using the AQCMD algorithm just described, we carried out preliminary QCMD simulations of liquid water and ice, using the q-TIP4P/F potential energy and dipole moment surfaces.qtip4pf We used sets of bond-angle centroids, defined as in Eq. (25) to describe the internal motion of the water monomers, and ‘Eckart-like’ rotational frames (see Sec. V.1) to describe the translations and overall rotations of the monomers. Numerical details and results of the simulations are presented in Sec. V.2.

V.1 Eckart-like frame

V.1.1 Properties of the frame

In molecular spectroscopy, the Eckart frame defines the orientation of a molecule by means of a reference geometry, usually taken to be the equilibrium geometry.eckart ; bright-wilson Here, we orient each water monomer using an ‘Eckart-like’ frame, in which the ‘molecule’ is the monomer ring-polymer, and the ‘reference geometry’ is the monomer quasi-centroid. This gives the conditions

[TABLE]

where we use the notation of Sec. III, such that $\mathbf{q}_{i}^{(\alpha)}$ are the atom polymer beads and $\overline{\mathbf{Q}}^{(\alpha)}$ are the atom quasi-centroids, with $\alpha=({\rm H}_{1},{\rm H}_{2},{\rm O})$. It is easy to show that these conditions are equivalent to

[TABLE]

where $\mathbf{Q}^{(\alpha)}$ are the atom centroids.

The first condition (Eq. (33a) or, equivalently, Eq. (34a)) thus places the centre-of-mass of the atom quasi-centroids at the centre-of-mass of the entire monomer ring-polymer, or equivalently, at the centre-of-mass of the atom centroids.

The second condition (Eq. (33b) or, equivalently, Eq. (34b)) specifies the orientation of the monomer ring-polymer about its centre-of-mass. To illustrate how this works, let us imagine that we have generated a new monomer quasi-centroid geometry $\overline{\mathbf{Q}}$ (using step A1 of the AQCMD algorithm), and that we now wish to define the ensemble of ring-polymers that ‘belong’ to this geometry. One can show that the second Eckart condition is equivalent to minimisingjorg ; kudin

[TABLE]

or, equivalently,

[TABLE]

In other words, each monomer ring-polymer in the ensemble is oriented so as to minimise its average mass-weighted distance from $\overline{\mathbf{Q}}$. In practice, this means that the atom centroids are usually very close to the atom quasi-centroids—see Fig. 6. The Eckart-like conditions thus ensure that the quasi-centroid-constrained ring-polymer distribution is compact with respect to the orientational degrees of freedom.minsp

V.1.2 Implementation

To implement the Eckart-like conditions, we carry out the propagation in cartesian coordinates $(\overline{\mathbf{P}},\overline{\mathbf{Q}})$, using the AQCMD algorithm of Sec. IV. This allows us to avoid using Euler angles or quaternions. Steps O2, B2, and A2 of the AQCMD algorithm are straightforward to implement by imposing the constraints of Eq. (25) and Eq. (34) on each monomer using SHAKE and RATTLE.ryckaerts ; andersen ; tuckerman Step B1 is implemented as follows.

First, to take into account that there are now external forces acting on the polymer beads, the expressions for the internal forces given in Eq. (30) are modified according to

[TABLE]

with

[TABLE]

where $\mathbf{q}_{i}^{(c)}$ is the centre-of-mass of the $i$-th replica, and $\mathbf{I}_{i}$ is its inertia tensor in the centre-of-mass frame. Having made this modification, we use Eqs. (29) and (30) to compute the internal forces on the quasi-centroid.

Second, we add the external forces on the monomer quasi-centroids, which are given by

[TABLE]

where

[TABLE]

with $\overline{\mathbf{Q}}^{(c)}$ the quasi-centroid centre-of-mass and $\mathbf{I}$ the quasi-centroid inertia tensor. The centre-of-mass force $\overline{\mathbf{f}}$ is easy to compute since it is just the average of the forces on the replica centres-of-mass

[TABLE]

but the torque $\overline{\bm{\tau}}$ is difficult to compute exactly.invert We therefore approximate $\overline{\bm{\tau}}$ as the average of the torques on the replicas

[TABLE]

This expression is consistent with the other approximations made to the forces (namely that $\overline{\bf P}$ is purely cartesian, and that the polymer springs can be neglected in the mean-field force), since Eqs. (39) and (41) agree to within second order in the displacements from the quasi-centroid. As with these other approximations, the reliability of Eq. (41) can be monitored by comparing QCMD static properties with the results of standard PIMD simulations.

V.2 Simulation results

V.2.1 Computational details

Condensed-phase water CMD, TRPMD and QCMD simulations were carried out in the three regimes studied in ref. 35, namely the compressed liquid at $600~{}$\mathrm{K}, the liquid at $300~{}$\mathrm{K}, and ice Ih at $150~{}$\mathrm{K}$$.

Unless noted otherwise, the CMD and TRPMD calculations used the same parameters as ref. 35, with periodic boundary conditionsfrenkel ; allen and simulation boxes of 128 molecules for the liquid and 96 for ice. At each temperature, a set of independent starting configurations was generated by propagating eight instances of the system classically for $100~{}$\mathrm{p}\mathrm{s} under the action of a global Langevin thermostat,[schneider-langevin ](#bib.bib75); [bussi-langevin ](#bib.bib76) followed by a $100~{}$\mathrm{p}\mathrm{s} TRPMD run under PILE-G (with the centroid relaxation time $\tau_{0}=100~{}$\mathrm{f}\mathrm{s}, and the non-centroid mode friction parameter $\lambda=0.5$).[pile ](#bib.bib65) The resulting eight equilibrated ring-polymer configurations were then propagated for a further $100~{}$\mathrm{p}\mathrm{s}, using PILE-G (with the same parameters as above) to give the TRPMD results, and the PA-CMD algorithm (with $\gamma=4$, $\tau_{0}=100~{}\mathrm{fs}$, and $\lambda=0.01$) to give the CMD results.

To compute the QCMD spectra for the liquid, we took the eight equilibrated configurations obtained after the first 100 ps PILE-G (see above), then propagated for $35~{}$\mathrm{p}\mathrm{s} using the AQCMD algorithm as described in Sec. [IV](#S4) and [V.1](#S5.SS1), with a global Langevin thermostat[bussi-langevin ](#bib.bib76) (time constant $\tau_{0}=100~{}$\mathrm{f}\mathrm{s}) acting on the quasi-centroids. The first $10~{}$\mathrm{p}\mathrm{s}$$ were discarded before calculating the dipole-derivative time-correlation functions.

For QCMD ice, we generated 15 independent equilibrated configurations (in the same manner as the 8 in the liquid simulations). These were then propagated with the AQCMD algorithm for $2.5~{}$\mathrm{p}\mathrm{s}, with a local Langevin thermostat[bussi-langevin ](#bib.bib76) acting on the quasi-centroids. For each of the 15 configurations, the system was then propagated under a global Langevin thermostat[bussi-langevin ](#bib.bib76) for three intervals of $2.5~{}$\mathrm{p}\mathrm{s}, separated by relaxation periods of $300~{}$\mathrm{f}\mathrm{s}, during which local thermostatting was used. The IR spectrum was then calculated by averaging over the three $2.5~{}$\mathrm{p}\mathrm{s} trajectories from the 15 independent simulations. Converged AQCMD results were obtained using $\gamma=8$ at 600 K, and $\gamma=32$ at 300 K. At 150 K, the calculation was stopped at $\gamma=128$, which has definitely converged the positions of the peaks and the $t\to 0$ limit of the dipole-derivative time-correlation function, but may not yet have converged the intensity of the stretch peak—see Appendix C. These values of $\gamma$ made the QCMD calculations approximately $2,\,8,\,\text{and}\ 32$ times more expensive than the corresponding CMD calculations, which used the PA algorithmhone2 ; mano-pacmd and thus smaller values of $\gamma$.

V.2.2 Overview of the QCMD spectrum

At 600 K, the QCMD spectra are in very close agreement with the TRPMD and CMD spectra (Fig. 7). Measured from the peak maximum, the QCMD OH stretch at $3600~{}\mathrm{cm}^{-1}$ is blue-shifted by about $15~{}\mathrm{cm}^{-1}$ relative to CMD and TRPMD. The bend at $1600~{}\mathrm{cm}^{-1}$ and the libration peak at $500~{}\mathrm{cm}^{-1}$ are almost identical in all three methods.

At 300 K, noticeable differences emerge between QCMD and the other two methods, due to the CMD curvature problem, and the broadening of TRPMD, although these problems are a lot less severe than in the gas-phase at this temperature.marx10 ; trpmd2 The CMD stretch peak is red-shifted relative to TRPMD by about $50~{}\mathrm{cm}^{-1}$; the QCMD stretch peak is blue-shifted by about $20~{}\mathrm{cm}^{-1}$. The bend peaks are in good agreement for all three methods. The libration band in QCMD is slightly red-shifted; the most possible cause for this small shift is the approximation made to the torques on the monomers in Eq. (41) (see Sec. V.2.4).

In $150~{}$\mathrm{K}$$ ice, there are major differences between the three sets of results, mainly in the stretch region. The curvature problem is strong, artificially broadening and red-shifting the CMD stretch by more than $100~{}\mathrm{cm}^{-1}$, and the TRPMD stretch is severely broadened. In contrast, the QCMD stretch is sharp, with resolved symmetric and antisymmetric bands. This resolution of the stretch bands is an artefact of the q-TIP4P/F surface and is also observed in classical simulations. Using a shorter time window of 500 $\mathrm{f}\mathrm{s}$ to coalesce the QCMD stretch peaks (not shown), gives a QCMD stretch that is blue-shifted by about $10~{}\mathrm{cm}^{-1}$ with respect to the TRPMD stretch, and is more than twice as intense. Note that (as mentioned above) the intensity of the QCMD stretch may not have converged with respect to $\gamma$. All three methods are in close agreement in the libration and bend regions, with a slight red-shift in the QCMD libration band which is smaller than the shift at 300 K (see Sec. V.2.4).

V.2.3 Stretch region of the QCMD infrared spectrum

To better assess the performance of QCMD in the stretch region of the spectrum ($>2000$ cm*-1*), we compare with two other methods: the coloured-noise version of TPRMD developed recently by Rossi et al.trpmd3 and the local monomer approximation (LMon) of Bowman and co-workers.lmon1 ; lmon2 ; trpmd2 Both methods are expected to do well in the stretch region of the spectrum, but to suffer from artefacts at lower frequencies.

In coloured-noise TRPMD, the white-noise PILE thermostat used in most TRPMD calculations (such as those reported above) is replaced by a generalised Langevin equation (GLE) thermostat that is designed to minimise the dynamical disturbance to the centroids at certain pre-tuned frequencies.trpmd3 Using the GLE(C) parametrisation of the thermostat in ref. trpmd3, (along with the other simulation parameters in this reference), we propagated eight independent 100 $\mathrm{p}\mathrm{s}$ TRPMD+GLE(C) trajectories at 300 and 150 K, using the i-PI package.ipi The resulting spectra in the stretch-region are shown in Fig. 8; (the lower-frequency parts of the spectrum, which are corrupted by the thermostat, are shown in Appendix D).

In the LMon method, the solid or liquid is equilibrated using standard PIMD techniques, then the Schrödinger equation for the nuclear dynamics is solved for each monomer independently, with all degrees of freedom except for the internal and few intramolecular modes of the monomer held fixed.lmon1 ; lmon2 ; lmon3 ; lmon4 Here, we make use of the results of previous LMon calculations,trpmd2 in which the PIMD equilibration was done using the same box size and simulation parameters as in the TRPMD calculations of Sec. V.2.1, and the internal modes of the water molecule were coupled to each of its three librational modes in turn, resulting in a four-mode approximation (LMon-4). We convolved the raw data obtained from these calculationsmariana with the same time-windows as used in the QCMD calculations (see Appendix A), to obtain the spectra shown in Fig. 8 for the stretch region at 300 and 150 K. The approximate treatment of the intermolecular modes means that LMon-4 gives a relatively poor description of the spectrum at 600 K (not shown) and of the libration at all temperatures (see Appendix D).

The agreement between the QCMD and TRPMD+GLE(C) stretch peaks is excellent: at 300 K the peaks are almost identical; at 150 K the TRPMD+GLE(C) peaks line up with the QCMD peaks, but are somewhat broader, with the bifurcation just visible; (the broadening is to be expected since the thermostat couples more strongly to the centroid at lower temperatures). There is probably some cancellation of errors between these two sets of results, as the QCMD method makes a number of static approximations (see Sec. II.2), and the GLE(C) thermostat interferes strongly with the dynamics of the librational modes which are coupled to the dynamics of the stretch modes. One piece of evidence for such errors is that the 300 K QCMD peak is slightly less intense than the TRPMD+GLE(C) peak, whereas one would expect the opposite to be true. Nonetheless the agreement between QCMD and TRPMD+GLE(C) in the stretch region is remarkable, suggesting that both methods give excellent approximations to the Matsubara-dynamics spectrum in this region.

The QCMD and LMon-4 results are also in good agreement in the stretch region (Fig. 8). Clearly one cannot use LMon-4 as a quantum benchmark, since it includes only a few degrees of freedom in the monomer dynamics (which is probably why it does not reproduce the bifurcation in the stretch peak). Nevertheless, we think the comparison in Fig. 8 highlights an important weaknesses in the QCMD approach (which is shared by CMDpaesnew and TRPMDtrpmd2 ). The LMon-4 spectrum gives a libration-bend combination band at roughly $2300~{}\mathrm{cm}^{-1}$ (see the Fig. 8 inset) and a Fermi resonance overlapping with the stretch peak at $3200~{}\mathrm{cm}^{-1}$; both these features are much more intense than the corresponding features in the QCMD, CMD and TRPMD spectra. While we cannot say how well the LMon-4 calculations have converged these features, it seems highly likely that QCMD, CMD and TRPMD grossly underestimate them. This is not surprising, as we expect such features to depend strongly on Matsubara fluctuations,mats ; mats-rpmd and also probably on real-time coherence—effects which QCMD, CMD and TRPMD omit.

V.2.4 Approximations to the quantum Boltzmann distribution

On the basis of the spectra reported above, it seems likely that QCMD gives a similarly good approximation to the quantum Boltzmann distribution for the condensed phase as it does for the two-dimensional model and the gas-phase. However, the condensed-phase QCMD calculations make extra approximations to the quantum Boltzmann distribution: Eq. (41) approximates the torques on the monomers, and the AQCMD algorithm only samples the exact quantum Boltzmann distribution in the limit $\gamma\to\infty$.

To test these static approximations, we compare (in Fig. 9) the AQCMD radial distribution functions (RDFs) for the liquid at 300 K and ice at 150 K with the exact RDFs (obtained from PIMD). As expected, the AQCMD RDFs are close to the exact values, although there are small errors at 300 K, indicating that AQCMD gives a slightly less structured liquid than it should. Given that the AQCMD RDFs are in almost perfect agreement with the exact values for 150 K ice, and that the intramolecular peaks in the O–H RDF (not shown) agree to graphical accuracy, we think that the (small) errors at 300 K are mainly the result of approximating the monomer torques by Eq. (41). There is perhaps also a very small contribution to the errors from incomplete convergence in $\gamma$; increasing $\gamma$ from 32 to 64 reduced slightly the errors in the RDFs, suggesting that higher values of $\gamma$ might reduce them further.

These results suggest that QCMD gives an excellent approximation to the quantum Boltzmann distributions in water and ice, except for small errors in the long-range correlations at 300 K, which are most likely due to the quasi-centroid torque estimator of Eq. (41). This estimator is also thought to be responsible for the small red shift in the libration band of the infrared spectrum at 300 K (see Fig. 7). Future work may thus be able to reduce both of these (small) errors, by developing better torque estimators.

VI Conclusions

By making minor approximations to the quantum Boltzmann distribution, we have obtained a method (QCMD) for propagating curvilinear mean-field centroids which is at least efficient enough to treat liquid water and ice on the q-TIP4P/F potential. For a two-dimensional model, the quantum Boltzmann approximations are found to be almost exact, and the infrared spectrum is close to that obtained from Matsubara dynamics, showing no sign of the ‘curvature problem’ that afflicts CMD at low temperatures. For liquid water (300 K) and ice (150 K), the errors made in the quantum Boltzmann approximations are found to be small, suggesting that the QCMD spectra are also very close to the Matsubara-dynamics spectra for these systems. The algorithm we developed to implement QCMD is slow (taking 32 times longer than PA-CMD for 150 K ice), but is do-able, and we hope to apply it soon to more realistic water potentials such as MB-pol, to allow comparison with experiment. It may be possible to speed up the algorithm by making it partially adiabatic, or perhaps to develop instead a curvilinear version of TRPMD.

As with other ring-polymer dynamics methods, such as CMD and (T)RPMD, QCMD is a (quasi-)centroid-following method, and thus cannot obtain information about the dynamics of the Matsubara fluctuations around the centroid. In general, QCMD will thus give poor results for non-linear operators waterpole (which depend explicitly on the dynamics of the Matsubara fluctuation modes), and for overtones and combination bands (which couple the centroid strongly to the fluctuation modes)—as an example of the latter, see Fig. 8.

How readily can QCMD be generalised to systems other than water? This question is equivalent to asking how readily collective bead coordinates can be identified, around which the ring-polymers are distributed compactly. Such coordinates are probably easy to find for liquids and solids made up of small, rigid molecules; but further work will be required to determine whether they can be found for systems containing floppy molecular units, such as solvated protons. Also, for reaction barriers below the instanton cross-over temperature, the collective bead coordinates will need to include the unstable mode of the instanton. A ‘holy grail’ would be to develop an algorithm that follows, automatically, the collective bead coordinate that minimises the spread of the ring-polymer distribution, thus giving the best mean-field approximation to Matsubara dynamics.

Acknowledgements.

We thank Mariana Rossi for sending us the LMon-4 data used to obtain Fig. 8. G.T. acknowledges a University of Cambridge Vice-Chancellor’s award and support from St. Catharine’s College, Cambridge. M.J.W. and S.C.A. acknowledge funding from the UK Science and Engineering Research Council.

Appendix A Details of the calculations of the infrared absorption spectra

The champagne-bottle potential used in Sec. II is

[TABLE]

where r=\big{(}q_{x}^{2}+q_{y}^{2}\big{)}^{1/2},\,r_{\mathrm{eq}}=1.8324, $D_{0}=0.18748$, and $\alpha=1.1605$ a.u. The mass of the particle is taken to be $m=1741.1$ a.u.

For the linear dipole-moment calculations (i.e. the 2D model and the condensed-phase water calculations), the infrared spectrum was obtained from the Kubo-transformed time-autocorrelation function (TAF) $C(t)$ of the cell dipole-moment derivative $\dot{\bm{\mu}}$tfmiller1 ; mano-pacmd

[TABLE]

where $V$ is the cell volume, $n(\omega)$ is the refractive index, $\sigma(\omega)$ is the absorption cross-section, and $f(t)$ is a window function that dampens the tail of the TAF, reducing the ringing artefacts.nr77 For the gas-phase water calculations (which employed the Partridge-Schwenke dipole-moment surface), we used the dipole TAF

[TABLE]

In two-dimensional and gas-phase calculations, where the cell volume $V$ is not defined, we take $\tilde{I}_{\dot{\bm{\mu}}}(\omega)$ or $\omega^{2}\tilde{I}_{{\bm{\mu}}}(\omega)$ to be “the spectrum” instead of $n(\omega)\sigma(\omega)$.

When applied to the exact quantum TAF, the damping envelope $f(t)$ serves the additional purpose of removing the recurrences that are caused by the real-time quantum coherence. In the two-dimensional model, the recurrences are well separated from the initial part of the TAF, allowing a sigmoid window with a sharp cut-off to be used,

[TABLE]

with parameters $t_{1/2}=400~{}$\mathrm{f}\mathrm{s}$,\,\text{and}\ \tau=25~{}$\mathrm{f}\mathrm{s}$$. In simulations of gaseous and condensed-phase water (Sec. III and V.2) the Hann windownr77

[TABLE]

is found to be more suitable. The cut-off time $\tau$ is set to 600 $\mathrm{f}\mathrm{s}$ for liquid water at 600 and 300 $\mathrm{K}$, and to 800 $\mathrm{f}\mathrm{s}$ for ice at 150 $\mathrm{K}$.

To obtain smooth LMon-4 spectra, for which the raw datamariana are a list of discrete transition frequencies $\{\omega_{i}\}$, we calculate

[TABLE]

where $\mu_{i}$ are the corresponding transition dipole moments and $\tilde{f}$ is the Fourier transform of the Hann window. The resulting spectrum is then scaled so that its area integrated between 2600 and 4500 $\mathrm{cm}^{-1}$agrees with QCMD.

Appendix B Treating the quasi-centroid momentum as purely cartesian

To obtain Eq. (15b), we assume that the quasi-centroid momentum $\overline{\bf P}$ contributes a cartesian term $e^{-\beta\overline{\bf P}^{2}/2m}$ to the mean-field-averaged partition function $Z_{p}(P_{R},L,R,\Theta)$, where $P_{R}$ and $L$ are the momenta conjugate to $R$ and $\Theta$. To investigate this approximation, it helps to specify $\Theta$. Let us use

[TABLE]

with $\theta_{i}=\arctan{y_{i}/x_{i}}$. We then obtain (without approximation)

[TABLE]

with

[TABLE]

To make $T(P_{R},L,\rho)$ factorise out of $Z_{p}(P_{R},L,R,\Theta)$ we need to approximate the moment of inertia $2m\rho^{2}$ by $2mR^{2}$, giving

[TABLE]

and to make the corresponding change $1/{\rho}^{\prime}\to 1/R^{\prime}$ in the integrand of Eq. (49). This gives

[TABLE]

Since

[TABLE]

where $\Delta r_{i}=r_{i}-R$, it follows that treating $\overline{\bf P}$ as purely cartesian is a good approximation if the ring-polymer distribution is sufficiently compact that we can neglect terms quadratic in ${\Delta r_{i}}$ in the moment of inertia in $T(P_{R},L,\rho)$.

Appendix C Convergence of AQCMD spectra for liquid water and ice

The convergence of the AQCMD infrared absorption spectra with respect to the adiabatic separation $\gamma$ at 600 and 300 K was tested by simulating liquid water using a box of 32 molecules, subject to periodic boundary conditions.frenkel ; allen The simulated spectra converged at $\gamma=8\ \text{and}\ 32$ respectively, as shown in Figure 10.

The convergence with respect to $\gamma$ for ice at 150 K was tested using a larger simulation box of 96 molecules (in order to give zero net dipole momenthayward ). The highest value tested was $\gamma=128$, which was sufficient to converge the positions of all the bands in the spectrum, but not the relative intensity of the bend and stretch bands (see Fig. 10). A value of $\gamma=64$ was sufficient to converge the dipole-derivative TAF at $t=0$. This, together with the abrupt nature of convergence at $300~{}$\mathrm{K}$$, suggests that the $\gamma=128$ spectrum is close to convergence at 150 K.

Appendix D LMon-4 and TRPMD+GLE spectra for water

For completeness we compare in Fig. 11 the full IR spectra calculated using the QCMD, TRPMD+GLE(C) and LMon-4 methods (the stretch regions of the spectra are also plotted in Fig. 8 and discussed in Sec. V.2.3). As mentioned in Sec. V.2.3, LMon-4 and TRPMD+GLE(C) are known to provide poor descriptions of the libration band.