Eva G Noya
email: eva.noya@iqfr.csic.es

Luis Sese

Rafael Ramírez

Carl Mcbride

Maria M Conde

Eva ; Noya

Eva G Noya

Luis M Sesé

Carlos Vega
email: cvega@quim.ucm.es

Path integral Monte Carlo simulations for rigid rotors and their application to water

Keywords: path integral simulations, water, thermodynamics, Monte Carlo

come   L'archive ouverte pluridisciplinaire

In 1948 Feynman proposed the path integral formulation of quantum mechanics 1,[START_REF] Feynman | Path-integrals and Quantum Mechanics[END_REF] . In the late seventies the works of Barker and of Chandler and Wolynes showed that this formulation could be implemented in the statistical mechanical study of condensed matter by performing classical simulations of a modified Hamiltonian [START_REF] Barker | [END_REF]4 . It was demonstrated that the partition function of a quantum system of N particles in its discretised form is formally identical to that of a classical system consisting of N ring polymers. Thus a number of the techniques and methods that had already been derived for classical simulations could now be adapted to perform quantum simulations. Since then path integral simulations have been used to study the behaviour of a large number of systems. A detailed description of the path integral technique in statistical mechanics, and its applications, can be found in several reviews [START_REF] Gillan | The path-integral simulation of quantum systems[END_REF]6 .

The path integral formulation can be implemented both in Monte Carlo (MC) and in molecular dynamics (MD) simulations. MC simulations only provide thermodynamic properties. Although there has been some controversy concerning the description of the correct dynamics, there are now MD methods, such as centroid molecular dynamics 7,8 or the ring polymer molecular dynamics 9 that approximate the quantum dynamics of the system.

The path integral method has also been extended to study rigid rotors. The first quantum MC simulations of rigid bodies were performed in the mid-eighties using a semi-classical approximation to derive the rotational contribution to the partition function 10,11 . This semi-classical approximation was also used within molecular dynamics in the beginning of the nineties 12 . A few years later Müser and Berne derived the exact propagator for the rotational contribution to the density function for rigid bodies 13,14 . Centroid molecular dynamics has also been recently extended to deal with rigid rotors 15 . These attempts to extend the path integral method to rigid bodies are motivated by the desire to describe the quantum behaviour of molecules in a computationally efficient way, since various condensed matter properties are more likely to be affected by inter-molecular vibrations rather than intramolecular contributions 14 . In this situation, it seems reasonable to describe the molecules as rigid rotors and ignore the intra-molecular vibrations. Usually intra-molecular vibrations exhibit high frequencies that require to use a large number of 'replicas' of the system, thus increasing considerably the computational cost of quantum simulations. This problem can be alleviated by the use of smart techniques such as the recently proposed ring polymer contractor method 16 . That said, it is our opinion that the description of molecules as rigid bodies is particularly interesting in its self; since it allows one to separate inter-molecular and intra-molecular quantum effects, thus assigning their relative influence on the different properties of the condensed phases.

In this work we shall describe the path integral formulation for rigid rotors proposed by Müser and Berne. Even though this formulation has already described in Refs. 13,14, the implementation of a path integral Monte Carlo involves many technical details which, in our opinion, are worth describing in detail. In addition, we show how this method can be applied to study quantum effects for a range of properties of water.

II. PATH INTEGRAL FOR RIGID ROTORS

A. General path integral formulation

The behaviour of a system of N quantum non-spherical rigid particles can be described by the Schrödinger equation:

ĤΨ λ = E λ Ψ λ (1)
where the Hamiltonian is given by:

Ĥ = T tra + T rot + Û (2)
where T tra is the translational kinetic energy operator of the centre of mass of the molecules, T rot is the orientational kinetic energy operator of the rigid molecules, and Û is the potential energy operator. The partition function for this system at inverse temperature β can be written as:

Z = λ Ψ λ |e -β Ĥ |Ψ λ = T r(e -β Ĥ ) = λ e -βE λ (3)
where Ψ λ and E λ are the eigenfunctions and the eigenvalues, respectively, of the Hamiltonian. The solutions of the Schrödinger equation form a complete orthonormal basis (i.e.

Ψ λ |Ψ λ ′ = δ λλ ′ ). The operator ρ = e -β Ĥ is defined by its Taylor expansion:

ρ = e -β Ĥ = ∞ k=0 Ĥk k! (4)
The second term in Eq. 3 implies first an integration over the coordinates of the system (r, ω)

followed by a sum over the states λ. Here r represents the set of Cartesian coordinates of the centre of mass of each particle and ω are the set of Euler angles that define the orientation of each particle in the system. Note that although we have used the eigenfunctions of the Hamiltonian to write the partition function, any complete basis set can be used due to the fact that the trace is invariant with respect to a change of basis [START_REF] Mcquarrie | Statistical Mechanics[END_REF] .

We now introduce the density function, which in the coordinate space is defined as:

ρ(r, ω, r ′ , ω ′ , β) = λ Ψ * λ (r, ω) exp(-βE λ )Ψ λ (r ′ , ω ′ ) (5)
The second term in Eq. 3 can be re-written as:

Z = T r(e -β Ĥ ) = λ Ψ * λ (r, ω) exp(-βE λ )Ψ λ (r, ω)drdω = ρ(r, ω, r, ω, β)drdω (6)
where one first sums over the states λ to obtain ρ(r, ω, r, ω, β) and then integrates over the coordinates. One interesting property of the density function is that the product of two density functions results in another density function:

ρ(r, ω, r ′′ , ω ′′ ; β 1 )ρ(r ′′ , ω ′′ , r ′ , ω ′ ; β 2 )dr ′′ dω ′′ = ρ(r, ω, r ′ , ω ′ ;

β 1 + β 2 ) (7)
Naturally this can be generalised to any given number of terms. Using this convolution property of the density function the partition function can be re-written as: Z = ρ r 1 , ω 1 , r 2 , ω 2 ; β P ...ρ r P -1 , ω P -1 , r P ω P ; β P ρ r P , ω P , r 1 , ω 1 ; β P dr 1 dω 1 ...dr P dω P

Usually it is not possible to solve the Schrödinger equation for a system of N rigid nonspherical interacting particles, so the eigenfunctions Ψ λ and eigenvalues E λ are unknown.

However, as mentioned before, the partition function (or what is the same the trace of the density matrix operator ρ) is invariant independently of the complete set of basis functions used. A convenient complete basis set is the one formed by the eigenfunctions of position and orientation operators which, for non-spherical particles, is given by:

|r t ω t = δ(r -r t )δ(ω -ω t ) (9)
Therefore, we can write the partition function using the eigenfunctions of the position and orientation operator by substituting the sum over the states λ by an integration over r and ω. In this way, the partition function can be calculated by evaluating the density function ρ t,t+1 (β/P ) ≡ ρ(r t , ω t , r t+1 , ω t+1 , β/P ) = r t ω t exp -β( T tra + T rot + Û )/P r t+1 ω t+1 (10) where |r t+1 ω t+1 are the eigenfunctions of the position operator. In principle, the exponential cannot be factorised because the potential and kinetic energy do not commute. However, using the Trotter formula [START_REF] Trotter | [END_REF] , one can write the following exact formula in the limit of infinite

P : lim P →∞ ρ(r t , ω t , r t+1 , ω t+1 , β/P ) = (11) r t ω t exp -β Û /2P exp -β( T tra + T rot )/P exp -β Û /2P r t+1 ω t+1
where P is the Trotter number. Given the ring polymer isomorphism, P is also known as as the number of replicas or beads.

The operator Û is diagonal in the coordinate representation, and for rigid rotors T tra and T rot commute. Therefore, for rigid rotors ρ t,t+1 (β/P ) can be approximated as:

ρ t,t+1 (β/P ) ≈ ρ t,t+1 pot (β/P )ρ t,t+1 tra (β/P )ρ t,t+1 rot (β/P ) (12)
where

ρ t,t+1 pot (β/P ) = exp - β 2P U (r t , ω t ) + U (r t+1 , ω t+1 ) , (13)
ρ t,t+1 tra (β/P ) = r t exp -β T tra /P r t+1 , (

) 14
ρ t,t+1 rot (β/P ) = ω t exp -β T rot /P ω t+1 . ( 15
)
where U (r t , ω t ) is the potential energy of N rotors whose positions and orientations are specified by r t , ω t .

It can be shown that the translational contribution is given by 5 :

ρ t,t+1 tra (β/P ) = M P 2πh 2 β 3/2 exp - N i=1 M P 2h 2 β (r t i -r t+1 i ) 2 , ( 16
)
where M is the rotor mass and r t i are the coordinates of the centre of mass of replica t of rotor i. In order to evaluate the partition function for a system of N free rotors it is necessary to evaluate the rotational propagator:

B. Rotational propagator for a free rotor

ρ t,t+1 rot (β/P ) = ω t exp - β P Trot ω t+1 (17)
which can be written as the product:

ω t exp - β P T rot ω t+1 = N i=1 ω t i exp - β P T rot i ω t+1 i ( 18
)
where |ω t i are the eigenfunctions of the orientation operator, which, as mentioned before, are Dirac delta functions:

|ω t i = δ(Ω -Ω t i ) (19)
Here Ω represents the three Euler angles (θ, φ, χ). To simplify the notation, we shall drop the subindex i, and from here on focus our attention on a single free rotor.

The eigenfunctions of the angular position |ω t can be expanded in a basis set of the eigenfunctions of the asymmetric top |JM K (the derivation of the eigenfunctions of the asymmetric top are given in Appendix A):

|ω t = J M K JM K|ω t |JM K (20)
Using this expansion, the rotational propagator can be rewritten as:

ω t exp - β P T rot J M K JM K | ω t+1 JM K (21)
As the functions |JM K are the eigenfunctions of the Schrödinger equation for the asymmetric top, it follows that:

exp - β P T rot JM K = exp - β P E (JM ) K JM K (22)
Using this, the rotational propagator can be written as:

ω t   J M K exp - β P E (JM ) K JM K | ω t+1   JM K (23)
Reordering this expression one has: The location of the laboratory frame defining the Euler angles is arbitrary. For convenience we choose a laboratory frame such that the Euler angles of replica t are all zero (i.e., Ω t = (0, 0, 0)). This change leads to:

ω t exp - β P T rot ω t+1 = J M K ω t | JM K exp - β P E (JM ) K JM K | ω t+1 (24)
|ω t = δ(Ω -Ω t ) = δ(Ω) (25)
and

|ω t+1 = δ( Ω -Ωt+1 ) (26)
A tilde is added to Ω (i.e., Ω) in order to remind ourselves that the Euler angles are defined in a laboratory frame in which the Euler angles of bead t are zero. Thus Ωt+1 are the Euler angles of bead t + 1 in this arbitrary frame.

To simplify this expression further, the eigenfunctions of the asymmetric top |JM K are expanded in a basis set formed by the eigenfunctions of the symmetric top (|JM K ):

|JM K = K A JM KK |JM K (27)
Eq. 24 can now be re-written as:

ω t exp - β P T rot ω t+1 = J M K δ(Ω) K A (JM ) KK |JM K dΩ exp - β P E (JM ) K K A (JM ) KK |JM K * δ( Ω -Ωt+1 )d Ω = J M K K A (JM ) KK Ψ JM K (0) exp - β P E (JM ) K K A (JM ) * KK Ψ * JM K ( Ωt+1 ) (28)
The eigenfunctions of the symmetric top are (see Ref. 19):

Ψ JM K (Ω) = 2J + 1 8π 2 1/2 exp(iM φ)d J M K (θ) exp(iKχ) (29)
where d J M K (θ) are the Wigner functions (given in Appendix B). Using this expression we have:

Ψ JM K (0) = 2J + 1 8π 2 1/2 d J M K (0) = 2J + 1 8π 2 1/2 δ M K (30)
and having made use of the following relation for Wigner functions: d J M K (0) = δ M K . Müser and Berne concluded that, substituting these expressions, Eq. 28 could be written as 13 :

Ψ * JM K ( Ω) = 2J + 1 8π 2 1/2 exp(-iM φt+1 )d J M K ( θt+1 ) exp(-iK χt+1 ) (31)
ω t exp - β P T rot ω t+1 = = J M K 2J + 1 8π 2 A (JM ) KM 2 exp - β P E (JM ) K d J M M ( θt+1 ) exp(-iM ( φt+1 + χt+1 )) ( 32
)
The propagator is a real quantity and this can be seen by symmetrising it with respect to M . This can be achieved by calculating the average of the M and -M contributions for any given J and K:

1 2 2J + 1 8π 2 A (JM ) KM 2 exp - β P E (JM ) K d J M M ( θt+1 ) exp(-iM ( φt+1 + χt+1 )) + 2J + 1 8π 2 A (J M ) K M 2 exp - β P E (J M ) K d J M M ( θt+1 ) exp(-i M ( φt+1 + χt+1 )) (33)
where M denotes -M . Given that (see Ref. 20):

d J M M (θ) = d J M M (θ) (34)
A (JM ) KM = A (J M ) K M (35)
and exp -

β P E (JM ) K = exp - β P E (J M ) K (36)
Eq.33 can be simplified to:

2J + 1 8π 2 A (JM ) KM 2 exp - β P E (JM ) K d J M M ( θt+1 ) cos(M ( φt+1 + χt+1 )) (37)
Using this result, we obtain the final expression for the orientational propagator:

ρ t,t+1 rot,i (β/P ) = ∞ J=0 J M =-J J K=-J f t,t+1 i,J,M, K exp - β P E (JM ) K (38)
where f t,t+1 i,J,M, K is a function of the relative Euler angles between beads t and t + 1 and is given by:

f t,t+1 i,J,M, K = 2J + 1 8π 2 A (JM ) KM 2 d J M M ( θt,t+1 i ) cos[M ( φt,t+1 i + χt,t+1 i )] (39)
C. Path integral simulations for a system of rigid rotors Once the three contributions to the density function (Eq.12), the potential (Eq. 13), the translational (Eq. 16) and rotational (Eq. 38) have been obtained, the partition function of a system of N rigid rotors can be calculated by substituting them into Eq. 8: If one uses a pairwise potential then U t = i j>i U (r t i , ω t i , r t j , ω t j ). One can now see that the partition function of a system of N quantum particles is isomorphic to that of a system of N classical ring polymers. Each replica t of molecule i interacts with the replicas with the same index t of the remaining particles through the inter-molecular potential U , and with replicas t -1 and t + 1 of the same molecule i through a harmonic potential, whose coupling parameter depends on the mass of the molecules (M ) and on the temperature (β = 1/k B T ), and through a term (ρ t,t+1 rot,i ) that incorporates the quantisation of the rotation and which depends on the relative orientation of replicas t and t + 1. The convention of the Euler angles used in the present calculations, as well as the conversion from Cartesian coordinates to Euler angles, are given in Appendix C. The procedure to obtain the Euler angles of replica t + 1 in the body frame of replica t is outlined in Appendix D.

Z = Q N V T = 1 N ! M P 2πβh 2
Note that the function ρ t,t+1 rot,i (β/P ) (Eqs. 38 and 39) depends solely on two angles, θt,t+1 and φt,t+1 + χt,t+1 . It is convenient to compute the density function for a grid of values over the angles θ and φ + χ and save these data in tabular form prior to any simulation. The value of the density function for any particular θ and φ + χ can then estimated using an interpolation algorithm in conjunction with the tabulated data.

The internal energy E can be calculated from the partition function Z using the thermodynamic relation:

E = - 1 Z ∂Z ∂β (41)
By substituting the partition function it can be shown that the internal energy can be calculated as:

E = K tra + K rot + U, (42)
where the functional forms of these three components are:

K tra = 3N P 2β - M P 2β 2 h2 N i=1 P t=1 (r t i -r t+1 i ) 2 , (43)
K rot = 1 P N i=1 P t=1 1 ρ t,t+1 rot,i    ∞ J=0 J M =-J J K=-J f t,t+1 i,J,M, K E (JM ) K e -β P E (JM ) K    , (44)
U = 1 P P t=1 U t , (45)
As with the rotational contribution to the density function, the numerator of the last term in Eq. 44 was calculated prior to simulation for a grid of angles θ and φ + χ and saved in tabular form.

The partition function for the N pT ensemble can be calculated using:

Q N pT = exp(-βpV )Q N V T dV (46)
The implementation of the N pT ensemble in PIMC has already been discussed in previous works [START_REF] Barrat | [END_REF]22 .

III. QUANTUM EFFECTS IN WATER

In this section we shall see how the path integral formulation for rigid rotors can be used to study quantum effects related to the atomic mass in water. Water exhibits such quantum effects even at room temperature. For example, properties such as the melting temperature, the temperature of maximum density (TMD), or the heat capacity at constant pressure, all show changes upon isotopic substitution [START_REF] Franks | Water: A Matrix of Life[END_REF][START_REF] Chaplin | [END_REF] . Experimental data for different water isotopes are given in Table I. In particular, the melting temperature of water is about 4.49K higher in tritiated water than in H 2 O. A more dramatic effect can be seen in the increase in the TMD, which is about 8.91K higher in tritiated water. In contrast, isotopic substitution of the oxygen mass has little influence on the properties of liquid water (see Table I). This indicates that quantum effects are mainly due to the light mass of hydrogen, which leads to small values of the principal moments of inertia. The lowest principal moment of inertia increases when hydrogen mass increases, whereas it remains almost unchanged when oxygen is substituted. The effect that this has on the magnitude of quantum effects can be seen by using the following approximate expression to estimate the quantum effects of a rigid asymmetric top [START_REF] Allen | Computer Simulation of Liquids[END_REF] : where F is the force that acts on the centre of mass of the particles, I A , I B and I C are the principal moments of inertia of the particles, and Γ A , Γ B and Γ C are the torques associated with each principal axis of inertia. This expression indicates that both the low moments of inertia as well as the strength of the hydrogen bond (which leads to high values of the average torque on the molecules) are responsible for the importance of quantum effects in water.

A -A cl N = h2 24(k B T ) 2 F 2 M + Γ 2 A I A + Γ 2 B I B + Γ 2 C I C - h2 24 cyclic 2 I A - I A I B I C (47)
The importance of nuclear quantum effects, as well as the limitations of performing classical simulations for water, have already been discussed by Stillinger and Rahman in 1974, one of the first articles published concerning simulations of water [START_REF] Stillinger | [END_REF] . In this pioneering paper it was pointed out that classical models are likely to overestimate the temperature difference between the TMD and the melting temperature of water, which experimentally is 3.98K. They speculated that this difference would be about 14K for classical models, which was obtained by plotting the difference between the TMD and the melting temperature for water, deuterated water and tritiated water, as a function of the inverse of the hydrogen mass and then linearly interpolating to infinite mass, which would correspond to the classical limit.

This hypothesis was also reiterated more recently by Guillot in his review of water models 27 .

A thorough investigation of classical rigid non-polarisable models eventually showed that the temperature difference between the TMD and the melting temperature was even larger; approximately 30K 28 .

The first simulations of water to explicitly include quantum effects were performed almost thirty years ago [10][11][12]29 . In these works, quantum effects in liquid water and deuterated water were investigated by using the rigid non-polarisable ST2 model. It was observed that liquid water is less structured than heavy water, and that classical water (which corresponds to the limit of infinite mass) is more structured than either of these isotopes. In the same period, Wallqvist and Berne studied quantum effects using a flexible model 30 .

Since these seminal works, a good number of publications have appeared treating the quantum effects in water, the majority focusing on the liquid phase 15,[31][32][33][34][35] . Some of these studies implemented potentials designed to be used within classical simulations. Even though this approach is valid when it comes to studying the effect that inclusion of quantum effects have on the diverse properties of water, a quantitative description requires the use of a potential that has been specifically designed to use within path integral simulations. So far a few water models have been developed expressly to be used within quantum simulations. These include both rigid 36 and flexible models that have either been fitted to experimental data 37,38 or to more accurate ab initio calculations 39 . Recently more sophisticated path integral Car-Parrinello molecular simulations of liquid water and ice I h have also been undertaken, which showed that ab initio calculations are also improved by incorporating nuclear quantum effects 40 . As far as we are aware, there are only a few studies that have considered quantum effects for ice I h 34, [41][42][43][44] . There are also a couple of studies that have examined the liquid-solid 35,38 and solid-vapour interface 35 .

In this work we treat the water molecule as being rigid. This means that our simulations are only capable of providing information about the low frequency inter-molecular librations (for water these are below 900cm -1 ), whereas high frequency intra-molecular vibrations (between 1500 cm -1 for bending and 3500 cm -1 for stretching) will be ignored. The number of replicas required to accurately reproduce the properties of a quantum system depends on the largest vibration frequency of the system:

P > hω max k B T (48)
From this expression it can be seen that at room temperature approximately P =30 replicas should be used for flexible models (for which the higher frequencies are around 3500 cm -1 ), whereas P =5 or 6 is sufficient for a rigid model (for which the higher frequencies are about 1000 cm -1 ). This permits a considerable reduction in the computational cost of the simulations.

F

IV. METHODOLOGY

In the quantum simulations presented in this work, water was described using the recently proposed TIP4PQ/2005 model 45,46 , which is the quantum counterpart of the TIP4P/2005 model 47 . The classical TIP4P/2005 model was found to provide the best overall description of water from among the many simple rigid non-polarisable models available in the literature [48][49][50][51] . In both models a Lennard-Jones (LJ) centre is located on the oxygen site, positive charges on the hydrogens and a negative charge along the bisector of the oxygenhydrogen vectors. The total energy of the system is given by:

U = i j>i 4ǫ σ r ij 12 - σ r ij 6 + m∈i n∈j q m q n r mn (49)
where r ij represents the distance between the oxygen atoms in molecules i and j, r mn is the distance between the charge q m of molecule i and charge q n of molecule j. σ and ǫ are the LJ parameters. The parameters of the TIP4P/2005 and TIP4PQ/2005 models are given in Table II. The only difference between the two models is that the hydrogen charges have been increased in the TIP4PQ/2005 model by 0.02e, which leads to stronger electrostatic interactions. This compensates for the loss of structure and the increase in energy that is observed when quantum simulations are performed. The same recipe was also used previously to obtain quantum counterparts of the classical SPC/F 52 and TIP5P 36 models. The increase in the charges enhances the dipole moment of the water molecule from 2.305D in TIP4P/2005 to 2.380D in TIP4PQ/2005 (higher multipole moments will obviously also change accordingly). Note that in some cases for flexible models the change in geometry caused by the incorporation of quantum effects with respect to the classical limit also leads to an enhancement of the dipole moment of the water molecule. Therefore, for flexible models it is generally not necessary to increase the charges in order to perform quantum simulations 38 .

In this work, the influence of quantum effects in water was investigated by perform- ing N pT PIMC simulations using the formulation for rigid rotors proposed by Müser and Berne 13,14 described previously. Classical N pT simulations were also performed for comparison. The simulation box contained 300 water molecules for the liquid phase, and 432 molecules for ices I h and II. The initial proton disordered configuration of ice I h was obtained using the algorithm of Buch et al. 53,54 . The LJ interaction was truncated at 8.5 Å. Standard long-range corrections for the LJ part of the potential were added. Coulombic interactions were calculated using Ewald summations. Simulations usually consisted of about 30,000 cycles for equilibration plus a further 100,000 cycles dedicated to obtaining averages. One MC cycle typically consisted of N P/2 Monte Carlo moves, N being the number of water molecules and P the number of replicas of the system. The configurational space was explored by using four types of movement attempts: translation of one single bead of one molecule (30%), rotation of a single bead of one molecule (30%), translation of a whole ring (20%) and rotation of all the beads of a given molecule (20%). The maximum displacement or rotation was adjusted in each case to obtain a 40% acceptance probability. After each N P/2 Monte Carlo moves one volume change was also attempted. The maximum volume change was adjusted to obtain a 30% acceptance probability.

V. RESULTS

Before presenting the results for bulk water, a preliminary check was undertaken to show that the rotational propagator was indeed able to reproduce the rotational kinetic energy of a free asymmetric top, for which the rotational energy can be analytically computed via evaluation of the partition function. The comparison between the simulation data with the analytical expansion is excellent both for water and tritiated water (see Fig. 1). It can also be observed that for temperatures above approximately 50K the quantum rotational energy of the free rotor is almost equal to the classical value (3/2k B T ), which means that for an isolated rigid water molecule quantum effects are only important below this temperature.

However, in condensed matter the situation is different because there are inter-molecular forces, in this case hydrogen bonds, that hinder the rotation of the molecules and lead to the appearance of noticeable quantum effects at much higher temperatures.

Before performing simulations for the liquid and solid phases, we need to choose the number of replicas, P , that are to be used in the simulations. The classical limit corresponds to one single replica, whereas the quantum limit is approached as the number of replicas tends to infinity. However, simulations can only be performed for a finite number of replicas.

In practice the number of replicas is chosen so that it is small enough that simulations do not become prohibitively expensive, but high enough so as to capture the main contribution of the quantum effects. P represents a compromise between statistical convergence and theoretical accuracy, so a study of how the desired property converges with the number of replicas has to be carried out. We examined the convergence of the potential energy and the total energy as a function of the number of replicas for liquid water at T =298K and p=1bar (see Fig. 2). It can be seen that both the potential and total energies increase with the number of replicas. There is a large increase for small number of replicas and then both magnitudes reach a plateau above 5 or 6 replicas. The quantum limit would be obtained by plotting the total energy as a function of the inverse of the number of replicas and taking the limit to infinite P . We found that the total energy is lower than the value at P → ∞ by about 3% for P =5 and by about 2% for P =7. In view of this we have chosen to use P =5 replicas at room temperature. Other authors have also used a similar number of replicas for water at room temperature 15,36,41,44 . For other temperatures, the number of replicas was chosen so as to keep the product P T approximately constant, i.e. taking P =5 at T =300K we arrive at P T ≈ 1500K.

A. Isotopic effects on the TMD and Cp

One of the idiosyncratic properties of water is the existence of a maximum in density. As mentioned previously, the location of the TMD is affected by variations in the hydrogen mass.

In particular, for deuterated water this maximum occurs 7K above the TMD of water and for tritiated water this increases to 9K. Therefore, we expect that nuclear quantum effects will shift the TMD of water to lower temperatures when quantum effects are implicitly incorporated. In addition, given that good water models reproduce the TMD of water For the system sizes studied in this work (N =300-360) the reaction field technique yielded slightly higher energies and densities than the Ewald summations, but the location of the TMD was unchanged 46 .

The results are shown in Fig. 3. The location of the TMD was obtained by fitting the simulation data to a quadratic or cubic polynomial. The TMDs obtained are given in Table III. The results show that, using quantum simulations, the TIP4PQ/2005 model predicts that the TMD of water occurs at 284(2)K, only 7K above the experimental result. We have seen that upon increasing the number of replicas the TMD shifts to 280(2)K 46 , indicating that simulation results become even closer to experimental results when more replicas are used. The results for deuterated water and tritiated water show that TIP4PQ/2005 is also able to qualitatively reproduce the shift to higher temperatures when the mass of the hydrogen isotope increases, in line with experimental observations. By considering the results for water, deuterated water and tritiated water (using P = 5 for the three molecules) one can see that the location of the TMD shifts 8K for deuterated water, and 12K for tritiated water, with respect to that of water. These are only slightly larger than the 7K and 9K, respectively, found experimentally. With regards to density, quantum effects affect differently the density depending on whether we are above or below the TMD. For high temperatures the number density increases with mass whereas for low temperatures the TMD the number density decreases when the mass is increased (see Fig. 3).

It is also interesting to study the shift in TMD when going from quantum water to classical water. Our simulations predict that the TMD can change as much as 30K when quantum effects are included. It seems that this is a typical shift for rigid non-polarisable models (a similar result was found for the TIP5P model 36 ). Similarly de la Pena, Razul and Kusalik estimated a shift in the melting point of ice I h for the rigid TIP4P model of about 35K when nuclear quantum effects are included 44 . However, the shift in the TMD between the quantum and classical limit could well be different for different types of potentials (for example, flexible or polarisable models). In particular, it has been found that for the flexible polarisable TTM2.1-F model 59 the explicit inclusion of quantum effects left the location of the TMD unchanged 60 . Further work is needed to clarify this. The TMD increases as the molecule becomes more and more classical (i.e. as the mass of the hydrogen isotope is increased). The number densities of water, deuterated water and tritiated water obtained from PI simulations along the room pressure isobar are shown in Fig. 3. Results from classical simulations along this isobar are also presented. As can be seen in Fig. 3 the number density at the maximum is hardly affected by the mass of the hydrogen isotope;

differences between the number densities at the maximum are within the estimated error bar. However, it seems that the number density at the maximum first decreases slightly on going from water to deuterated and then tritiated water, and then increases a little in the classical limit. Experimentally it has been observed that the number density at the maximum decreases by ≈ 0.7% on going from water to tritiated water (see Table I). PIMC simulations, to a lesser extent, also reflect this decrease (≈ 0.2%). We have shown in previous work that the internal energy of water is non-linear when plotted as a function of the inverse of the mass of the hydrogen isotope 46 . The same also seems to be true here for the density of water at the maximum. This indicates that the behaviour of classical water cannot be obtained from a simple extrapolation of results obtained for water, deuterated water and tritiated water.

Another interesting calculation is the difference between the TMD and the melting temperature obtained via classical and quantum simulations (∆T = T T M D -T melt ). It has been found that for classical simulations of rigid non-polarisable models the TMD is situated about 30K above the melting temperature (i.e., ∆T =30K) 28 , which is much higher than the ∆T =4K found experimentally. Quantum simulations for the flexible q-TIP4P/F model predict that the difference between the TMD and the melting point is also about 30K 38 .

Preliminary direct coexistence simulations 61,62 of the melting temperature together with the TMD calculations presented before indicate that the difference might be about 20-22K for the TIP4PQ/2005 model, which improves upon the classical prediction, but that is still far from the 4K found experimentally. This suggests that although the inclusion of nuclear quantum effects reduces the value of ∆T other features of real water need to be incorporated in the model, such as polarisability, in order to quantitatively reproduce the experimental difference.

From the simulations performed along the 1 bar isobar it is straightforward to evaluate the heat capacity at constant pressure for water as well as other water isotopes (C p = ∂H ∂T p , H being the enthalpy). It has been found experimentally that the heat capacity of liquid water is considerably affected by the isotopic substitution of the hydrogen atom (see Table I). In particular, the heat capacity is about a 10% higher for deuterated water than for water at room temperature 63 water using the same model 64 , which again indicates that C p is significantly affected by the inclusion of quantum effects. One might think that a good description of the heat capacity of liquid water at room temperature could be obtained by using a classical description with a model in which quantum effects are implicit though the parametrisation of the model.

However, to the best of our knowledge, none of the rigid non-polarisable models proposed so far has been able to provide a quantitative description of the heat capacity of liquid water 64 .

As an example, results for TIP4P/2005 are also shown in Fig. 4, which are in poor agreement both with experiments and with quantum simulations. This suggests that quantum effects need to be incorporated explicitly in order to obtain a quantitative description of the heat capacity of water. This is further supported by the finding that PIMC simulations with TIP4PQ/2005 also reproduce the experimental heat capacity of ice I h from very low temperatures up to room temperature (results for ice I h are shown later). These results demonstrate that the main contribution can be captured using a rigid model and that the contribution from the intra-molecular degrees of freedom is small. This is not unexpected; intra-molecular vibrations exhibit very high frequencies (≈ 3000 cm -1 ) so at room temperature only the ground state is populated and therefore there is little or no contribution to C p from the intra-molecular vibrations.

B. Equation of state of ices

A significant deficiency of classical simulations is the inability to reproduce the equation of state of solids at low temperatures. One of the consequences of the third law of thermodynamics is that the thermal expansion coefficient should tend to zero at zero temperature 67 , which is equivalent to saying that the density should remain constant at low temperatures.

However, in previous work, we have seen that classical simulations of both TIP4P/2005 47 and TIP4P/Ice 68 are unable to reproduce the correct curvature of the equation of state of ices at low temperatures 69 . This was attributed to the fact that quantum effects become increasingly important as the temperature decreases. To see whether the description of ice at low temperatures could be improved when quantum contributions were incorporated, we are shown in Fig. 5. These results show that a good agreement with the experimental data can be obtained when quantum contributions are explicitly included; the equation of state has the same curvature as the experimental curve, now in concordance with the third law of thermodynamics. The same was also found to be true for ice I h 45 Both classical and quantum simulations were performed for ice II at T =100K. The results are given in Fig. 8. The first peak of the oxygen-oxygen distribution function is considerably lower for quantum simulations with TIP4PQ/2005 than for classical simulations with TIP4P/2005. As far as we know, as yet there are no experimental data for the atomic distribution function of ice II at this thermodynamic state, but it is expected that quantum simulations provide a better description of the structure at low temperatures (as was the case for ice I h 45 , the only ice for which the oxygen-oxygen atomic distribution function has been experimentally measured). The effects are much larger if one examines the oxygen-hydrogen and hydrogen-hydrogen distribution functions, for which differences between classical and quantum simulations extend further than the first peak. In particular, the classical hydrogen-hydrogen distribution function exhibits a large number of well defined peaks which become considerably smoother in the quantum limit. The results for ice I h and II also suggest that nuclear quantum effects might affect more significantly the hydrogen-hydrogen atomic distribution function in proton ordered ices, such as ice II.

D. Thermodynamic coefficients for ice I h

In previous work, we demonstrated that classical simulations using both the TIP4P/2005 and the TIP4P/Ice models were unable to provide a good description for many thermodynamic coefficients for ice I h 69 . In particular, it was shown that classical simulations resulted in a poor description of the heat capacity at constant pressure and of the thermal expansion coefficient. The thermal compressibility, on the other hand, was described reasonably well.

We checked whether the description of some of these thermodynamic coefficients could be improved by performing quantum PIMC simulations with the TIP4PQ/2005 model. The heat capacity at constant pressure and the thermal expansion coefficient can be calculated from the simulations that trace out the room pressure isobar (data shown in Fig. 5). The heat capacity was obtained by fitting the enthalpy to the function H = a + bT 2 + cT 3 and differentiating this fit with respect to the temperature. The isothermal compressibility was calculated by performing simulations at p =-500, -250, 0, 250 and 500 bars for temperatures between 100K and 250K. The density along each of these isotherms could be nicely fitted to a straight line. The isothermal compressibility was computed by differentiating the density with respect to the pressure from the fit (κ T = 1 ρ ∂ρ ∂p T

). Once the thermal expansion coefficient and the isothermal compressibility are known, the pressure coefficient (β V = ∂p ∂T V ) can be readily computed via β V = α/κ T . The thermal coefficients as obtained from PIMC simulations using the TIP4PQ/2005 model, as well as those form classical MC simulations using the TIP4P/2005 model 69 are shown in Fig. 9. The thermal coefficients derived from the experimental equation of state of Feistel and Wagner are also shown for comparison 65 . The results indicate that, except for the thermal compressibility, for which a similar accuracy is obtained in quantum and classical simulations, quantum simulations considerably improve the description of the thermodynamic coefficients. In particular, the heat capacity at constant pressure, which was not reproduced by classical simulations at any temperature, is now nicely reproduced from room temperature all the way down to zero kelvin. This is in line with the results for liquid water, where again the experimental heat capacity could only be reproduced by the explicit inclusion of quantum effects. These results strongly suggest that quantum effects are crucial Finally, the description of the pressure coefficient (β V ) is also considerably improved when including quantum effects. This simply reflects that the thermal expansion coefficient is improved in a quantum description of the system, since β V = α/κ T . The isothermal compressibility is little affected by quantum effects; almost all the change in the pressure coefficient is due to a good description of the thermal expansion coefficient.

F

E. Relative energies of ices at zero kelvin

Finally we have also computed the relative energies between various ice phases at zero temperature. It has been found that some classical water models result in a rather good de- scription of the phase diagram of water 47,[START_REF] Abascal | [END_REF]75 . However, there is still room for improvement.

For example, it has been found that usually ice II is over-stabilised with respect to ice I h , for some models, so much so that ice II completely removes ice I h from the phase diagram 76 .

A preliminary outline of the phase diagram for a particular model can be obtained by estimating the coexistence pressures between the competing solid phases at zero temperature.

At zero kelvin phase transitions occur with zero enthalpy change, so a calculation of entropy is avoided. Assuming that the change in energy and density between two solid phases is almost independent of pressure at zero kelvin (which is indeed a rather good approximation for ices), the calculation of coexistence pressures between two ices at zero temperature can be estimated from 77 :

p = - ∆E ∆V p=0 (50)
Therefore, by simply calculating the energy and density of the solid phases at zero temperature and zero pressure one can obtain a reasonable estimate of the coexistence pressure at zero temperature.

The properties at zero temperature were computed for ices I h , II, III, V and VI. Empty hydrates structures sI, sII and sH, which have been shown to be the stable solid phases at negative pressures 78 , were also considered. Simulations were performed along the zero bar isobar in the temperature range from 250K to 77-100K. The energy at zero temperature was obtained by fitting the data to the function E = a + bT 2 + cT 3 , from which one can estimate E(T = 0K). The energies obtained using this procedure are represented in Fig. 10. Energies are given relative to the energy of ice I h , which experimentally is the most stable phase at zero temperature and at zero pressure. The results show that both classical MC simulations using TIP4P/2005 and quantum PIMC simulations using TIP4PQ/2005 predict that ice I h is the most stable phase, in agreement with experimental results. It is also observed that the relative energies of ices II, III, V and VI change when quantum effects are explicitly taken into account. In particular, ice I h is destabilised with respect to ice II by about 0.2kcal/mol, so that now the relative stability of ice II with respect to ice I h is much closer to the experimental value. This indicates that quantum effects are larger in ice I h than in ice II, resulting in a de-stabilisation of the former. Ices III, V, and VI are also stabilised with respect to ice I h , but to a smaller extent (by about 0.1 kcal/mol). Finally, the relative stability of the empty hydrate structures are not appreciably changed. Taking everything into account, we can identify three different families of ices according to the importance of quantum effects. The first family includes ice I h and the empty hydrate structures sI, sII and sII, which are influenced the most by quantum effects. The second family will be that formed for ices III, V, and VI, and finally, ice II, which is the least affected by quantum effects, forms the third family.

The reason why quantum effects make distinct contributions to the various ice phases can be understood by looking at the geometrical arrangement of the four molecules that form a hydrogen bond with a central one. These molecules form a nearly perfect tetrahedron in ice I h and a slightly deformed tetrahedron with deviations of about 10 degrees in the hydrates.

On the other hand the deviations from the perfect tetrahedron are of about 30 degrees for the ices II, III, V and VI. As a result, the strength of the hydrogen bond is larger in ice I h and the hydrates than in the remaining ices and, as can be seen in Eq. 47, this results in an increased value for the average forces and torques, boosting the quantum influences.

Once the energies at zero temperature have been calculated, the coexistence pressures can also be obtained. The results are given in Table IV. As can be seen, the coexistence pressure between ices I h and II is the most affected coexistence line. It decreases from about 2090 bar in classical simulations (TIP4P/2005) to about 195 bar in quantum simulations (TIP4PQ/2005), resulting in a much better agreement with the experimental results. The coexistence lines between the remaining ices (II, III, V and VI) are also affected, but to a lesser extent. In general, quantum results are closer to the experimental data than the results obtained using classical simulations. Finally, coexistence pressures between I h and the hydrates or between hydrates themselves are largely unaffected by quantum effects.

In summary, phase transitions between solid phases belonging to the different families are the most affected by the inclusion of quantum effects, whereas phase transitions involving phases of the same family are in general less affected. Therefore, the explicit inclusion of quantum effects is crucial if one wishes to reproduce the phase transitions between solid phases belonging to different ice families, especially the transition I h -II, whereas usually phase transitions between ice phases of the same family can be calculated by means of classical simulations in conjunction with a good classical model.

VI. CONCLUSIONS

In this work it has been shown in detail how the formulation of the path integral for rigid rotors, derived by Müser and Berne 13 , can be applied to water. Using this formulation, a large number of properties of liquid water, ices and hydrates have been studied by PIMC simulations using the TIP4PQ/2005 model, which was specifically designed to use within quantum simulations 45 . For liquid water, isotopic effects on the TMD and on the heat capacity have also been considered. The results show that in general a better description of water is obtained when quantum effects are included, although some properties can also be reasonably described in classical simulations with a good classical model. In addition, PIMC simulations with TIP4PQ/2005 reproduce the experimental isotopic effects on the TMD and the heat capacity of liquid water. Quantum effects have been found to be crucial when it comes to reproducing many properties of water and ices. In particular, quantum effects have been found to be most important with regards to the properties of ices at low temperatures, which is not entirely unexpected.

In particular, classical simulations fail to reproduce the curvature of the equation of state at low temperatures found experimentally and imposed by the third law of thermodynamics 69 .

Our results show that this can be corrected by including quantum effects resulting in physically agreeable equations of state for ices I h 45 and II over a quite broad range of temperatures.

Classical simulations overestimate the first peak in the oxygen-oxygen radial distribution function of ice I h at 77K, which again is brought into to agreement with experiment when quantum simulations are performed 45 . The same behaviour was found for ice II and is expected to occur for other ices. As a result of the better description of the equation of state the coefficient of thermal expansion at low temperatures of ice I h is also improved when quantum effects are included. In addition to the properties of ices at low temperatures, quantum effects have also been found to be important when it comes to reproducing the heat capacity of ice I h and water at all temperatures. In addition, isotopic effects on the heat capacity of liquid water have also been captured.

We also found that the magnitude of quantum effects is different for different ices and, therefore, they need to be included if one wishes to improve the description of phase transitions. In particular, it has been found that ices can be classified into three different families, according to the importance of quantum effects: the first family is formed by ice I h and the hydrates structures sI, sII and sH, for which quantum effects are the largest, the second family comprises ices III, V and VI, and the third family is formed by ice II, for which quantum effects are the smallest. Phase transitions between ices belonging to different families change when quantum effects are included, whereas transitions between ices belonging to the same family are only slightly affected by quantum effects. As quantum effects are also different for liquid water and ice I h , the melting point of water is also affected by quantum effects; it has been found in previous works that the melting point shifts to lower temperatures when quantum effects are included 38,44 . In general, quantum affects should also affect any property that involves two phases for which quantum effects are different. For example, in previous work it has been shown that quantum effects improve the description of the enthalpy of vaporisation 46 and the sublimation enthalpy 71 .

However, it has been found that other room temperature properties, although also af- fected by quantum effects, can be properly described using a classical model. This can be explained because classical models are usually fitted to reproduce some experimental data at room temperature, so in some way quantum effects at this temperature are implicit in the model. For example, it has been found that the structure of liquid water and ices I h and II above 250K is reproduced with similar accuracy in classical and quantum simulations (although only quantum simulations can reproduce isotopic effects) 46 . The densities of ices can also be reproduced at room temperature with good accuracy using classical simulations 45 .

The isothermal compressibility of ice I h shows little improvement with respect to classical simulations at all temperatures, which also might be related to the fact that quantum effects are influenced little by pressure.

With regards to the TMD, several classical models were proposed that reproduce the experimental TMD 47,55 , although the temperature difference between the TMD and the melting temperature is largely overestimated (they usually predict a 30K difference 28 instead of the 4K found experimentally). Preliminary PIMC calculations indicate that this difference might be reduced for TIP4PQ/2005 to ≈ 22K, which means that, although some improvement is achieved, other features of real water need to be included, such as polarisability and flexibility, in order to obtain a better agreement with experiment.

By using a rigid model, we have ignored the influence of quantum effects in the intramolecular degrees of freedom. Despite this seemingly drastic approximation the results presented here seem to indicate that for many properties the main quantum contributions arise from the inter-molecular degrees of freedom. Comparison between quantum simulations of rigid and flexible models will be be very useful to quantify the relative importance of quantum effects on inter-molecular and intra-molecular degrees of freedom. Therefore, the eigenfunctions of the Hamiltonian will also be eigenfunctions of L2 and Lz , but not of Lc :

L2 Ψ = J(J + 1)h 2 Ψ J = 0, 1, ..., ∞ (54)
Lz Ψ = M hΨ M = -J, ..., 0, ..., J

The solutions for the Schrödinger equation for an asymmetric top will be denoted as |JM K .

The integer K is not a true quantum number (i.e., it does not quantise any observable) it is simply a number used to label the (2J + 1) possible values of the energy available for each value of J and M . The functions |JM K can be obtained by expanding them in a basis set formed by the eigenfunctions of the spherical top (|JM K ):

|JM K = K A JM KK |JM K (56)
Using the bca convention, the energies E ) can be obtained solving the following secular determinant (see Refs. 79,80) for each value of J and M :

H KK = 1 2 (B + C)[J(J + 1) -K 2 ] + AK 2 H KK±2 = 1 4 (B -C)[J(J + 1) -K(K ± 1)] 1/2 [J(J + 1) -(K ± 1)(K ± 2)] 1/2 (57)
where K ranges from -J to +J. The remaining elements of the determinant are zero. A, B and C are the rotational constants, A = h 4πIa , B = h 4πI b and C = h 4πIc . Note that since I a ≤ I b ≤ I c , it follows that A ≥ B ≥ C. This determinant has dimensions of (2 × J + 1). Therefore, (2 × J + 1) eigenvalues are obtained for each value of J and M , which are, in general, all different. These (2×J +1) energy levels are labelled with the K index. However, as M does not appear in the determinant, there is a (2 × J + 1) degeneracy in the energy associated with M . Note that if we had chosen a different convention for the reference system, e.g. the abc, the secular equation that we would need to solve would be different

     X Y Z      =     
cos φ cos θ cos χsin φ sin χcos φ cos θ sin χsin φ cos χ cos φ sin θ sin φ cos θ cos χ + cos φ sin χsin φ cos θ sin χ + cos φ cos χ sin φ sin θ

-sin θ cos χ sin θ sin χ cos θ           x b y b z b      (59)
Let us assign a particular body frame to the molecule of water. Taking the bca convention, the body frame of the water molecule is chosen so that the b principle axis (i.e., that associated with I b ) is assigned to x, c is assigned y and a is assigned to z (remember that

I a ≤ I b ≤ I c ).
Under this convention, the b axis lies along the H-O-H bisector, the c axis is perpendicular to the molecular plane and the a axis is parallel to the line connecting the two hydrogen atoms. With this choice the coordinates of the oxygen are (α, 0, 0) (the centre of mass in water is located upon the H-O-H bisector, slightly below the oxygen atom, and we shall define the x direction such that α is positive). The coordinates of the hydrogen atoms are be given by (β, 0, -γ) and (β, 0, γ). Notice that the value of β should be negative, since α was taken to be positive. We denote

(X O , Y O , Z O ), (X H 1 , Y H 1 , Z H 1 ) and (X H 2 , Y H 2 , Z H 2 )
as the coordinates of the oxygen and of the hydrogens, respectively, in the laboratory frame, whose origin is fixed to be the centre of mass of the molecule. As mentioned before, the centre of mass of the molecule coincides with the origin of the body frame fixed in the molecule. For the oxygen, the coordinates in the body frame fixed in the molecule (α, 0, 0) are related to those of the laboratory frame through Eq. 59. The multiplication of matrices in Eq. 59 leads to three equations:

X O = α(cos φ cos θ cos χ -sin φ sin χ) (60)
Y O = α(sin φ cos θ cos χ + cos φ sin χ)

Z O = -α sin θ cos χ (61)
Analogously, the coordinates of the hydrogens in the laboratory frame are given by:

X H 1 = β (cos φ cos θ cos χ -sin φ sin χ) -γ cos φ sin θ (63)
Y H 1 = β (sin φ cos θ cos χ + cos φ sin χ) -γ sin φ sin θ (64)

Z H 1 = -β sin θ cos χ -γ cos θ (65)
X H 2 = β (cos φ cos θ cos χsin φ sin χ) + γ cos φ sin θ

Z H 2 = -β sin θ cos χ + γ cos θ (67)
In summary, we have obtained 9 equations having 5 unknowns (φ and χ go from 0 to 2π and, therefore, their value can only be unambiguously obtained from the knowledge of both their sine and cosine; whereas θ varies from 0 to π so that it is only necessary to know its cosine). The Euler angles for an instantaneous configuration where the atomic coordinates are (X 0 , Y 0 , Z 0 ), (X H

cos θ = Z H2 -Z H1 2γ (69)
Subtracting Eqs.66 from 63:

cos φ = (X H 2 -X H 1 ) 2γ sin θ (70)
and subtracting Eqs. 67 from 64:

sin φ = (Y H 2 -Y H 1 ) 2γ sin θ (71)
Finally, the Euler angle χ can be obtained from Eq.62:

cos χ = -Z O α sin θ (72)
and adding Eq.63 and Eq.66:

sin χ = cos φ cos θ cos χ - (X H 1 + X H 2 ) 2β 1 sin φ (73)
or, alternatively, adding Eq.64 and Eq. 67:

sin χ = -sin φ cos θ cos χ + (Y H 1 + Y H 2 ) 2β 1 cos φ (74)
In the special case that θ = 0 or θ = π the expressions given above are not valid because the denominator vanishes, resulting in a singularity. The probability of obtaining exactly θ = 0 or θ = π is very small during a simulation. In these special cases, alternative expressions can be obtained by evaluating of the rotation matrix for the particular value of θ. When θ = 0, sin θ = 0 and cos θ = 1 and, therefore, the rotation matrix (Eq. 59) becomes i.e, in this case, the rotation can be seen as a simple rotation about the z-axis by an angle φ ′ = φ + χ. In the case θ = 0 there is no a unique way of defining φ and χ individually, as any combination of φ and χ having the same value of φ ′ (φ ′ = φ + χ) will provide the same final configuration. Here for the particular case θ = 0 we decided to assign φ = 0 , and with this choice the sine and cosine of the χ angle can be readily obtained using a procedure similar to that outlined above for a general case but using the rotation matrix given in Eq.

     cos(φ + χ) -sin(φ + χ) 0 sin(φ + χ) cos(φ + χ) 0 0 0 1      (

75.

Using this procedure we obtain :

sin(χ) = Y H 1 β (76)
and

cos(χ) = X H 1 β (77)
Finally, when θ = π, sin θ = 0 and cos θ = -1 so that the rotation matrix is now:

     -cos(χ -φ) sin(χ -φ) 0 sin(χ -φ) cos(χ -φ) 0 0 0 -1      (78)
In this case, this rotation can be seen as a simple rotation about the z-axis by an angle φ ′ = χ -φ. Again it is not possible to assign in a unique way values of χ and φ. For this reason we arbitrarily assigned in this case φ = 0 so that χ is obtained as:

sin χ = Y H 1 β (79) cos χ = - X H 1 β (80)
Appendix D. Obtaining the relative Euler angles of replica t + 1 with respect to those of replica t.

Let us focus on two replicas t and t + 1 of a certain molecule. Suppose that R t are the laboratory frame instantaneous coordinates of a certain site of replica t with respect to its centre of mass (the centre of mass of replica t). These coordinates can be obtained from the body frame coordinates of that site (r b ) using the rotation matrix M t : Inverting the previous equation one obtains

R t = M t r b (81)
r b = M t -1 R t (82)
The relative Euler angles of replica t + 1 with respect to those of replica t can be computed by expressing the instantaneous coordinates of replica t + 1 (R t+1 ) in the reference system of replica t by:

R ′ t+1 = M t -1 R t+1 (83)
In other words, replica t+1 is rotated using the rotation matrix M t -1 to obtain its orientation with respect to that of replica t. The relative orientation of molecule t + 1 with respect to that of molecule t is given by the atomic coordinates R ′ t+1 . The Euler angles associated with this orientation, i.e., the Euler angles of replica t + 1 with respect to those of replica t, can be computed using the procedure described in Appendix C. Physics in the special issue: "EIGHTH LIBLICE CONFERENCE on the statistical mechanics of liquids". We would like to thank the referee for their useful suggestions.

The answer to referee's criticism and the changes made to the manuscript are detailed below.

Answer to the referee:

1) Page 13, line 37: The text here refers to an "increase in the dipole moment". Since water does NOT have the structure of a dipolar fluid, it would seem that it is indeed increases in the higher moments that lead to the behaviour discussed here. The text needs to be appropriately corrected. (One way would be in use language of 'increased electrostatic interactions'.)

We have rewriten this sentence. We now clarify that the change in the charges changes not only the dipole moment but also the higher order multipoles.

2) Page 13, line 59: Need to clarify that it is the LJ potential that is truncated at 8.5 Å.

We have clarified that is the LJ potential that is truncated at 8.5 Å.

3) Page 17: Sentence starting on line 12, appears to cast some doubt on the results/conclusions here (e.g. line 20-23). This statement needs to be clarified -i.e. a 4

K improvement relative to what? what is the sensitivity to isotope?

We have rewritten this to clarify that the TMD of quantum water changes a little bit when going for P=5 to P=7. The agreement with respect to experimental results using P=7 is quite good. This is now clearly stated on the paper. if the authors wish to discuss these results.

We have now specified that the estimate of the melting temperature is obtained by direct coexistence simulations.

5) Page 18, Fig. 3:

There is an inconsistency in the density data in this figure, namely the values of the densities do not appear to be converging to the classical results as the hydrogen mass is being increased. Is there a simple explanation here? (e.g. perhaps these simulations were otherwise conducted differently?) This is an important point that needs clarification.

We now state clearly in the paper that we are plotting number densities rather than mass densities. The differences in density at the maximum are within the error bar for water, deuterated water, tritiated water and classical water. In previous work we found that the internal energy of water treated classically can not be obtained by extrapolating the results of water , deuterated water and tritiated water. From the results of this work it is clear that the same is true for the densities along the room pressure isobar. The reason for the failure of the extrapolation is that the properties of water do not change linearly with the inverse of the mass of the hydrogen. This point is now introduced and discused in the paper. We agree with the referee that Ref. 42 is appropriate in that part of the text and accordingly we have added the citation.

We hope that we have properly addressed all the comments raised by the referee.

Best wishes

Dr Eva G Noya

FIG. 1 :

 1 FIG. 1: Kinetic rotational energy of the isolated H 2 O (filled red squares) and 3 H 2 O (open blue circles) molecule as a function of temperature. There is a good agreement between the path integral simulations and the rotational energy obtained from the theoretical partition function of an asymmetric top for H 2 O (solid line) and 3 H 2 O (dashed line) geometry. For clarity, the rotational energy of 3 H 2 O has been shifted 0.1kcal/mol in the y-axis.

FIG. 2 :

 2 FIG. 2: Convergence of the potential energy U (circles) and the total energy E (squares) as a function of the number of replicas P for liquid TIP4PQ/2005 water at T=298K and p =1bar. The dotted and dashed lines are only guides to the eye.

FIG. 3 :

 3 FIG. 3: Isotopic effects on the TMD of water along the room pressure isobar. Number densities (i.e number of molecules per unit of volume) as a function of temperature at room pressure are presented.

.

  The heat capacities obtained from PIMC simulations with TIP4PQ/2005 model for liquid water, deuterated water, tritiated water as well as classical water simulated with the TIP4P/2005 model are shown in Fig. 4. It can be seen that quantum simulations with the TIP4PQ/2005 model are able to quantitatively reproduce the heat capacity of liquid water for a broad range of temperatures. In addition, the increase in the heat capacity upon isotopic substitution is also quite accurately captured. The agreement between simulations and experimental data for deuterated water is quite remarkable. Classical simulations with the TIP4PQ/2005 model are about 25% higher than those for quantum

  performed PIMC simulations using the TIP4PQ/2005 model in order to obtain the equation of state of ice II. The results of quantum simulations with TIP4PQ/2005 together with the results of classical simulations with TIP4P/2005 and the experimental data of Fortes et al.70

FIG. 4 :

 4 FIG.4: Heat capacity at constant pressure for water and other water isotopes as a function of temperature at p=1bar calculated by means of PIMC simulations with the TIP4PQ/2005. The heat capacity for classical water was simulated with the TIP4P/2005. Experimental results for water65 (filled circles) and deuterated water 66 (filled squares) are also shown for comparison.

FIG. 5 :

 5 FIG. 5: Equation of state of ice II at p =1bar as calculated from PIMC simulations with TIP4PQ/2005 and from classical MC simulations with TIP4P/2005. Experimental data are also shown for comparison 70 .

FIG. 6 :

 6 FIG. 6: Atomic distribution function of liquid water at 298K and 1 bar as calculated from classical MC simulations with TIP4P/2005 (blue dashed line) and quantum PIMC simulations with TIP4PQ/2005 (red solid line). Experimental data (black dotted line) are also shown 74 .

FIG. 7 :

 7 FIG. 7: Atomic distribution function of ice I h at 250K and 1 bar as calculated from classical MC simulations with TIP4P/2005 (blue dashed line) and quantum PIMC simulations with TIP4PQ/2005 (red solid line)

FIG. 8 :

 8 FIG. 8: Atomic distribution function of ice II at 100K and 250K and at 1 bar as calculated from classical MC simulations with TIP4P/2005 (blue dashed line) and quantum PIMC simulations with TIP4PQ/2005 (red solid line)

FIG. 9 :

 9 FIG. 9: Thermodynamic coefficients (κ T , α, β V and C p ) of ice I h at p=1bar, as calculated from classical simulations with TIP4P/2005 and from quantum simulations with TIP4PQ/2005. Experimental data are also shown for comparison 65 .

Page 27 of 44 URLFIG. 10 :

 4410 FIG. 10: Estimate of the relative energies at zero temperature of the ice phases as calculated from classical MC simulations with TIP4P/2005 and quantum simulations with TIP4PQ/2005 model. These data are taken from Refs. 45,71,76,78. Experimental data taken from Ref. 77 is also given for comparison.

FIG. 11 :

 11 FIG.11:Reference system fixed in the water molecule adopted in this work. This convention is usually denoted as the bca convention, because the x,y and z axis coincide with the b, c and a principal axes of inertia (thus I a < I b < I c ). The origin of the coordinate system is located at the centre of mass of the water molecule.

  In an asymmetric top all three moments of inertia are distinct (I a = I b = I c ). In this situation the Hamiltonian commutes with L2 and with Lz , but not with Lc :[ L2 , Ĥrot ] = 0 [ Lz , Ĥrot ] = 0 [ Lc , Ĥrot ] = 0 (53)

(

  see, e.g., Refs. 79-81).

Y H 2 =

 2 β (sin φ cos θ cos χ + cos φ sin χ) + γ sin φ sin θ

  to submit a revised version of the manuscript entitled "Path integral Monte Carlo simulations for rigid rotors and their application to water" by E. G. Noya, L. M. Sese, R. Ramirez, C. McBride, M. M. Conde, and C. Vega, for publication in Molecular

4 ) 1 Page 43 of 44 URL

 4144 Page 17, line 57. The nature of these "preliminary simulations" needs to be clarified

6 )

 6 Page 30, line 9. The work in Ref. 42 (see section III.C) concluded that "This difference of about 35 K should represent an upper bound to the expected shift of the thermodynamic melting temperature." It would seem a citation to this estimate is appropriate here.

TABLE I :

 I Experimental data for different water isotopes. The heat capacity at constant pressure is given at T=290K. H 2 O 2 H 2 O 3 H 2 O 1 H 2 18 O

	T melt (K)	273.15 276.97 277.64 273.43
	TMD (K)	277.13 284.33 286.55 277.35
	ρ T M D × 10 2 (molec./ Å3 ) 3.344 3.326 3.322 3.347
	TMD-T melt (K)	3.98 7.36 8.91	3.92
	C p (cal mol -1 K -1 )	17.93 19.60	-	-

1

TABLE II :

 II Parameters of the models TIP4P/2005 and TIP4PQ/2005.

	TIP4P/2005 3.1589 93.2	104.52	0.9572 0.1546 0.5564
	TIP4PQ/2005 3.1589 93.2	104.52	0.9572 0.1546 0.5764

Model

σ( Å) ǫ/k B (K) HOH (deg) d OH d OM q H (e)

  tions was reduced by using the reaction field method56 , rather than Ewald summations, to account for the long range electrostatic forces. It has been shown that Ewald summation and reaction field provide very similar results for liquid water57 and for Stockmayer fluids58 .

	Molecular Physics
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	47,55 , it would be interesting to check whether the TIP4PQ/2005 model is also able to reproduce O n l
	the experimental TMD.	y
	The equations of state of water, deuterated water, tritiated water and classical water
	were calculated at p =1bar. N pT PIMC simulations were performed at six different temper-
	atures in each case. As the maximum in density is reflected in the third significant figure of
	the density, especially long simulations are required to reduce the statistical error, so each

simulation consisted of at least 3 million MC cycles. The computational cost of the simula-Page 16 of 44 URL: http://mc.manuscriptcentral.com/tandf/tmph

TABLE III :

 III Temperature of maximum density at p =1bar as obtained from PIMC simulations of the TIP4PQ/2005 model for several water isotopes.

			System	TIP4PQ/2005 Expt.
		H 2 O (P=5)		284(2)	277.13
		2 H 2 O (P=5)		292(2)	284.33
		3 H 2 O (P=5)		296(2)	286.55
		Classical H 2 O		307(2)	-
	F r o	3.40
	)
	P 3 3.35 2 (molec./ Å e 3.30 ρ 10 e H 2 O Classical H 2 O 2 H 2 O r 3 H 2 O
		240 3.25	R 260 280	T (K) 300	320	340
			e
				v i e
					w
						O n l
						y

  and for hydrate sI71 .It is usually found that classical simulations using simple models of water tend to overestimate the height of the first peak in the oxygen-oxygen distribution function for both liquid water47 and for ice I h69 . It is well known that quantum effects lead to less structured liquids15,29,72,73 and solids41 and so one might think that this can be corrected by including quantum effects. Radial distribution functions for liquid water and ices I h and II obtained

	Page 21 of 44	Molecular Physics
	C. Structure of ices
	from classical simulations for TIP4P/2005 and quantum simulations with TIP4PQ/2005 at
	relatively high temperature (i.e. T >250K) and room pressure are presented in Figs. 6, 7 F and 8. Differences are visible although relatively small, being larger for the oxygen-hydrogen o r and hydrogen-hydrogen distribution functions, which is not unexpected since quantum ef-
	P
	e
	e r
		R
		e
		v i e
		w
		O n l
		y
		21

URL: http://mc.manuscriptcentral.com/tandf/tmph fects are mainly due to the hydrogen mass. The situation is different at low temperatures.

Page 28 of 44 URL: http://mc.manuscriptcentral.com/tandf/tmph

		Molecular Physics
	Phase TIP4P/2005 TIP4PQ/2005	Expt.
	I h -II	2090	195	140(200)
	I h -III	3630	2727	2400(100)
	II-V	11230	15731	18500(4000)
	II-VI	8530	10935	10500(1000)
	III-V	3060	1998	3000(100)
	V-VI	6210	6848	6200(200)
	I h -sI I h -sII o F r I h -sH	-4174 -3379 -4072	-3948 -3249 -3933	---
	sII-sI P e 2787 sII-sH -7775	2267 -7557	--
		e r
		R
		e
		v i e
			w
				O n l
				y

TABLE IV :

 IV Coexistence pressures (in bar) at zero temperature obtained from quantum PIMC simulations with TIP4PQ/2005 and classical MC simulations with TIP4P/2005. These data are taken from Refs. 45,71,76,78. Experimental data taken from Ref. 77 is also given for comparison.

Page 31 of 44 URL: http://mc.manuscriptcentral.com/tandf/tmph

  I. N. Levine, Molecular spectroscopy (John Wiley & Sons Inc., New York, 1975). 82 C. G. Gray and K. E. Gubbins, Theory of molecular fluids (Clarendon Press, Oxford, 1984).

	Page 35 of 44	Molecular Physics Molecular Physics Molecular Physics	Page 32 of 44
	80 P. R. Bunker and P. Jensen, Molecular symmetry and spectroscopy (NRC Research Press, Ot-
	tawa, 1998).
		(+y)
	F o r F o r F o r F o r F o r
	P P P P P
	e e e e e
	e r e r e r e r e r
		R R R R R
		e e e e e
		v i e v i e v i e v i e v i e
		w w w w w
		O n l O n l O n l O n l O n l
		y y y y y
		35

URL: http://mc.manuscriptcentral.com/tandf/tmph 81 URL: http://mc.manuscriptcentral.com/tandf/tmph

  1 , Y H 1 , Z H 1 ) and (X H 2 , Y H 2 , Z H 2 ) can therefore be obtained by solving the set of equations defined by Eqs. 60 -68. Subtracting Eq. 68 from 65 we obtain:

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgments

This work was funded by grants FIS2010-16159, FIS2010-15502 and FIS2009-12721-C04-04 of Dirección General de Investigación, S-0505/ESP/0299 of Comunidad Autónoma de Madrid and 200980I099 of CSIC. Useful discussions with Prof. J.L.F. Abascal are gratefully acknowledged.

Appendix A. Asymmetric top eigenfunctions

The rotation of an asymmetric free rotor can be obtained by solving the Schrödinger equation:

Ĥrot Ψ = T rot Ψ = EΨ (51) where T rot is the Hamiltonian associated with the angular momentum:

In this equation Lx , Ly and Lz are the three components of the angular momentum and I xx , I yy and I zz are the three components of the momentum of inertia. To solve this equation it is convenient to choose a reference system so that the x, y and z axis are located along the three principal axes of inertia, denoted as a, b and c. We adopt the convention that

Note that there is no unique way to identify x, y, z with a, b, c (see Refs.

79,80). For example, one could associate x with a, y with b, and z with c, which is usually referred to as abc convention. Alternatively we could choose to identify x with b, y with c and z with a, which is usually known as bca convention, and it is this convention that was used in this work (see Fig. 11). The choice of axis is highly relevant since it defines the Euler angles that appear in the three components of the angular momentum.

Appendix B. Wigner functions

The Wigner functions are given by 82 :

where the sum over χ is restricted to those values that do not lead to negative factorials.

Appendix C. Obtaining the Euler angles of a water molecule from the

Cartesian coordinates of its interaction sites

In our implementation of the MC algorithm we used Cartesian coordinates. This is convenient when it comes to computing the potential energy between two water molecules.

However, to evaluate the density function of the quantum rotational energy, the orientation of the molecules must be specified in terms of Euler angles. Therefore, we need a procedure to obtain the Euler angles that define the orientation of a given molecule from the Cartesian coordinates of its interaction sites. We used Euler angles (θ, φ, χ) as defined in Ref.81; θ varies from 0 to π and φ and χ go from 0 to 2π. Let us denote X, Y, Z as the orthogonal axes of a laboratory frame fixed in the space. Let us assume that we associated three orthogonal axes to the molecule, namely x, y, z, which define the body frame (with its origin located at the centre of mass of the molecule). We shall assume that both sets of orthogonal axes are right handed. The orientation of the molecule can be defined by three Euler angles.

The three Euler angles are defined by the operations required to move the molecule from an initial orientation, where x, y, z are coincident with X, Y, X (having both set of axes a common origin) to its current configuration (where the two set of axes also have a common origin). Rotations are counterclockwise. First, a rotation is performed about the Z-axis by an angle φ, so that the axes X, Y, Z change to x ′ , y ′ , Z. Secondly, we perform a rotation by an angle θ about the y ′ -axis obtained from the previous rotation (x ′ , y ′ , Z → x ′′ , y ′ , z ′ ).

Finally, a new rotation is performed around the z ′ -axis of the frame fixed in the body (

x ′′ , y ′ , z ′ to x, y, z). Using this convention, the coordinates in laboratory frame of a site of the molecule, X, Y, Z (i.e., R ) can be obtained easily from the coordinates of that site in body frame, x b ,y b ,z b (i.e., r b ) using the expression :