A Stereographic Projection Path Integral Study of the Coupling between the Orientation and the Bending Degrees of Freedom of Water

A Monte Carlo path integral method to study the coupling between the rotation and bending degrees of freedom for water is developed. It is demonstrated that soft internal degrees of freedom that are not stretching in nature can be mapped with stereographic projection coordinates. For water, the bending coordinate is orthogonal to the stereographic projection coordinates used to map its orientation. Methods are developed to compute the classical and quantum Jacobian terms so that the proper inﬁnitely stiff spring constant limit is recovered in the classical limit, and so that the nonconstant nature of the Riemann Cartan curvature scalar is properly accounted in the quantum simulations. The theory is used to investigate the effects of the geometric coupling between the bending and the rotating degrees of freedom for the water monomer in an external ﬁeld in the 250 to 500 K range. We detect no evidence of geometric coupling between the bending degree of freedom and the orientations. © 2008 American Institute of Physics . (cid:3) DOI: 10.1063/1.2925681 (cid:4)


I. INTRODUCTION
The path integral 1 approach to statistical mechanics has become a tool of choice for the investigation of quantum effects in clusters and other types of condensed matter at finite temperatures.  Deste a number of recent advances, the majority of path integral simulations have focused on atomic systems with strict adherence to Cartesian coordinates.This limitation is technical in nature, as the simple remapping of a Euclidean space by curvilinear coordinates greatly complicates both the formal and the numerical aspect of Monte Carlo path integral ͑MCPI͒ methods.However, strict adherence to Cartesian coordinates makes the simulation of molecular condensed matter a formidable task.In a recent article, 24 we use a relatively simple harmonic model for condensed matter and analytical, finite Trotter number ͑k m ͒, solutions of the path integral to investigate its convergence properties as a function of the mode frequencies and temperature.We find that when values of the spring constants typical of covalent modes alternate with values of the spring constants typical of intermolecular modes, there exist temperature ranges inside which the convergence of the path integral is highly nonuniform even at relatively elevated values of k m .We find that for a difference in the values of the spring constants as small as two orders of magnitude, constraining the high frequency degrees of freedom ͑DF͒ increases the numerical efficiency of the path integral substantially.This pattern is observed with both linear as well as cubically convergent solutions.
8][29] In fact, the Feynman quantization is generally possible, though more chal-lenging, in differential manifolds even when spaces have configuration dependent curvature 27 and torsion. 29These complications have been of little concern to molecular physicists since most of the non-Euclidean spaces generated by holonomic constrains have constant curvature and no torsion.1][32][33][34][35][36][37][38] There remain several numerical difficulties for the development and implementation of MCPI methods in curved spaces such as the imposition of spacial periodic boundary conditions on the time evolution propagator and the presence of multiply connected spaces.These difficulties have forced investigators to use relatively expensive alternatives to random walks like the vector-space MCPI method [31][32][33][34][35][36] or the use of fixed axes approximations. 30It is difficult to define the path integral measure when using open sets and periodic boundaries and special methods have to be developed to evaluate the path integral even in the most trivial cases. 29In a series of recent articles 24,[39][40][41][42][43] we have introduced and tested methods based on stereographic projection coordinates ͑SPCs͒ that can overcome the difficulties associated with the boundary conditions.The SPC map ⌽ : R d ← M d is a bijection from the manifold to an equidimensional Euclidean space where the coordinates range from negative to positive infinity.In contrast, mapping of typical d-dimensional manifold M d is achieved with open sets, as, for example, ellipsoids of inertia are mapped with traditional Euler angles or with quaternions.Open set maps do not cover the entire manifold.For example, the angular variable cannot include all the points in a circle of radius R; the point at = 0 and 2 must be excluded since maps must be single valued.The common remedy is to patch the map at a point ͑or set of points with zero Riemann measure͒ by using periodic boundary conditions.With a single SPC map there remains one point that is not covered in the manifold, but it is at infinity.
Using SPC maps, it is also possible to expand paths using random series, after generalizing the Feynman-Kaç formula in manifolds.This latter advance allows us to develop fast converging algorithms without the need to evaluate the gradient or the Hessian of the potential, and to employ efficient estimators based on numerical derivatives of the action. 20,21The numerical algorithms are based on the DeWitt formula, 27 but make use of a one-to-one SPC map between points in the d-dimensional manifold M d and an equivalent Euclidean space.The extension of these techniques to non-Euclidean manifolds mapped with SPCs has given new insights into the quantum effects of two important hydrogen bonded systems. 24,43ur recent work 24 with water clusters provides the next challenge for the developments of SPC-MCPI theory of condensed molecular matter.Constraining all the internal DF of the water molecules produces a very efficient and accurate approach for the simulation of water clusters in the 100Ͻ T Ͻ 250 K range.However, above 250 K the bending mode of water may contribute to the heat capacity, while below 500 K the hydrogen-oxygen stretching modes are predominantly in the ground state.Therefore, there is a clear need to extend our SPC-MCPI formalism to allow for simulations of molecules that include relatively soft internal DF.
The purpose of the work we report in the present article is to develop the necessary differential geometry using coordinates amenable for path integration in a challenging case, the bending water molecule with rigid stretches.It is our intention to use the formalism to study carefully the nature and the amount of coupling that takes place between the bending and the other DF.The particular space for the bending translating and rotating water molecule requires the development of the classical and quantum Jacobian terms so that the proper infinitely stiff spring constant limit is recovered in the classical limit, and so that the nonconstant nature of the Riemann Cartan curvature scalar is properly accounted in the quantum simulations.
For the classical Jacobians, the difference between the infinitely stiff spring limit Euclidean space simulation and the curved subspace simulation arising from holonomic constraints has been explored in detail using as an example the bending degree of freedom of a triatomic specie. 44The Lagrangian and the Hamiltonian for these two physical models are identical; however, there is a difference in the classical Boltzmann distribution.The difference is the result of expressing the partition function first in the flat space followed by the transformation of variables and the application of the infinite force constant͑s͒ limits, as opposed to the alternative procedure of expressing the Boltzmann distribution after the transformation and the application of the holonomic constraints.Frenkel and Smit 44 have argued that the former procedure is the correct one.In the curved manifolds one evaluates the metric tensor in the ͑3n − c͒-dimensional space ͑n atoms and c constraints͒, whereas in the c infinitely stiff springs model one remaps the R 3n space with internal DF and then evaluates the Jacobian from the metric tensor in the 3n flat space.Let us denote the metric tensor obtained by curvilinear remapping of R 3n spaces as g ͑E͒ , and the generic metric tensor in a non-Euclidean space as g .Let g 1/2 and g ͑E͒ 1/2 represent the square root of the determinant of the two metric tensors, i.e., the Jacobians; then, the volume element dV ͑E͒ , obtained by remapping R 3n and taking the infinite spring constant limit is in general different than the same in the curved manifold.Moreover, the ratio g ͑E͒ 1/2 / g 1/2 may still depend on the 3n − c coordinates, yielding possibly different results for the simulated properties with these two models.
For the quantum term of the Jacobian, the holonomic transformation of the action is straightforward and can be carried out with the usual tensor analysis machinery.However, as it has been pointed out by Kleinert, 29 the path integral measure in non-Euclidean spaces has to be derived by slicing in the flat space first, and then applying the transformation of coordinates to the resulting discretized measure.In a space with curvature, the order in which slicing and coordinate remapping is performed produces different expressions for the path integral measure.Kleinert produces a new quantum equivalence principal, which determines how a path integral in flat space behaves under holonomic ͑and even non-holonomic͒ transformations.The end result of Kleinert's work is a unique expression for the quantum Jacobian term, usually interpreted as a "quantum potential" that, in the absence of torsion, takes the following form

R, ͑2͒
where R is the Riemann-Cartan curvature scalar.For the majority of applications to molecular dynamics the curvature is a constant term and it can be ignored.For the bending trimer with constrained stretching, however, the curvature is a function of configuration and must be included in path integral simulations.This article is structured as follows.In Sec.II the formal development of the metric tensors associated with our choices of rotations and bending coordinates is presented in detail.In Sec.II we also develop in some detail the expressions for the classical and quantum Jacobian.The results of our numerical tests designed to study the couplings between the orientation and the bending of water in an external field are presented in Sec.III.Section IV contains our conclusions.

II. METHODS
In this section we derive two equivalent formulations for the metric tensor of a translating, rotating, and bending nonlinear triatomic, with fixed stretching coordinates that are suitable for path integral simulations.SPCs are developed to map both the rotation and the bending DF.We discover orthogonality properties when using Euler angles for the rotation, Cartesian coordinates for the translations, and a SPC for the bending degree of freedom.The bending coordinate is orthogonal to the rotation and to the translation ones.The same orthogonality properties remain when we remap the ellipsoid of inertia submanifold with SPCs.These properties are not essential for the simulations with path integrals but do provide the opportunity to decrease computer time costs when simulating molecular clusters.Additionally, orthogonality properties simplify the formulation of theories derived by splitting the propagation operator.Such theories will be the subject of future investigations.In the course of the analysis necessary to develop the metric tensors we find relatively simple analytical proofs for the following important property of the g ͑E͒ 1/2 / g 1/2 ratio: It is generally true that shown that the introduction of reference body-fixed configurations allows one to treat separately the rotations and translations from the internal DF.

A. The rigid three-body problem and the body-fixed frame
Consider a series of maps with the associated Jacobian matrices, metric tensors, and kinetic energy expressions that are useful for the cases when the intramolecular DF are characterized by infinitely stiff spring constants.We find that the best approach to arrive at a set of coordinates amenable for path integration is to develop the mapping into two stages.The first stage is to transform from the Euclidean R 3n space mapped by Cartesian coordinates to the Euler angles.
For a molecule with n atoms the map ⌽ 1 is defined by the transformation where , , are the three Euler angles defined by the element R of the rotation group, R = cos cos − cos sin sin sin cos + cos cos sin sin sin − cos sin − cos sin cos − sin sin + cos cos cos cos sin sin sin − cos sin cos ,

͑4͒
and r i is a vector in R 3 associated with atom i.The c subscript refers to the center of mass vector.R can be used in two ways.One can either rotate the axes of the body frame ͑BF͒, a rotating non-inertial frame in which the inertia tensor is diagonal, or rotate the reference body-fixed configuration.
When we express the metric tensor associated with ⌽ 1 −1 we drop the subscript "E" since ⌽ maps points of a Euclidean R 3n space into a six-dimensional non-Euclidean ͑and generally curved͒ subspace.Specifically, for the single water molecule example, ⌽ 1 is a map from a nine-dimensional Euclidean space mapped with Cartesian coordinates to a Cartesian product of two subspaces ⌽ 1 : R 9 → R 3 I 3 .There are several possible ways of defining ⌽ 1 , depending on a multitude of choices available for the definition of the body frame axis, and the choice between passive versus active rotations of the body-fixed axis.For the water molecule we choose the fol-lowing three vectors to represent the reference body frame axis, in units of bohr.Two items are important at this point.First, the configuration in Eq. ͑5͒ is defined so that the molecule is in the center of mass frame, i.e., and so that the inertia tensor is diagonal,

͑7͒
To transform from the body-fixed to the laboratory frame one rotates this rigid configuration and translates the center of mass.⌽ 1 −1 can be represented with a single set of equations, where R is the 3 ϫ 3 matrix in Eq. ͑4͒ and Once the map is specified and verified to be one-to-one, the task is to transform the metric tensor g from its form in Euclidean spaces mapped with Cartesian coordinates, to its form within the manifold.If we agree to represent raised indices as the column labels associated with Cartesian coordinate x , and the lower ones as row labels associated with the generalized coordinate q , then the Jacobian matrix ‫ץ‬x / ‫ץ‬q = J associated with ⌽ −1 can be represented by a 6 ϫ 3n matrix.One can use the transformation of tensors law which, with the row and column convention for the Jacobian matrix that we introduce here, can be expressed as a matrix product G ' =JGJ T .The metric tensor in Eq. ͑9͒ contains the effective masses associated with the coordinates.In the case of Euclidean spaces mapped by Cartesian coordinates, the effective mass is the physical mass of the body located by the coordinate point; in other spaces ͑or in Euclidean spaces mapped with curvilinear coordinates͒ the effective mass may become configuration dependent.With this convenient definition one can write a general formula for the kinetic energy K in generalized coordinates ͑with no velocity dependent interactions͒, To find expressions for the element of the Jacobian matrix let us introduce some additional notation.Let ‫ץ‬ R represent a set of 3 ϫ 3 matrices containing the derivative of the elements of the rotation matrix R with respect to q .There are a number of useful properties for the set ‫ץ‬ R, some of which are immediately obvious, This simple result yields trivial elements of the Jacobian matrix of ⌽ 1 −1 ,

͑13͒
In the last expression we reserve the Roman index i to label the atom number.The Roman indices k , kЈ , kЉ are used in what follows to span the R 3 subspace associated with r i .The second expression for the Jacobian matrix elements is With our notation we can obtain a generic expression for the transformation of an n body metric tensor like Eq. ͑9͒, where the fact that the Cartesian n body metric tensor is diagonal has been used.We must inspect three cases separately: where we exchange the order of summation and we use the fact that the body-fixed frame is in the center of mass.This is an important and general result; simply stated, it is the reflection that the translational and rotational coordinates as expressed by the map in Eq. ͑8͒ are always orthogonal.

͑18͒
The symbol ⌫ represents a set of nine tensors that operate in the R 3 subspace associated with the body-fixed Cartesian coordinates for atom i, The properties of this set are of general importance and we shall investigate them further.Since the set is symmetric, i.e., ⌫ = ͓⌫ ͔ T , as one can verify, we only need the expression for six matrices,

͑25͒
The transformation of Eq. ͑9͒ using Eq.͑5͒ for the expressions of r i ͑BF͒ , and Eqs.͑16͒ and ͑18͒, yields the following metric tensor: Using Eq. ͑7͒, and after some trivial manipulations, one can show that this metric tensor can be written as follows, In the Appendix we demonstrate that the metric tensor takes the following general form: The metric tensor in Eq. ͑28͒ can be obtained by exchanging with and then permuting the rows and columns corresponding to the d and d base vectors in Eq. ͑27͒.We have used both the passive and the active form in past computations. 24However, both expressions contain the same physics, and the difference between the passive and the active map and the resulting kinetic energy expression is immaterial.In Sec.II B, we transform Eq. ͑27͒ into its equivalent with SPCs, rather than Eq.͑28͒.However, from the nature of the SPC map it becomes evident that the difference between the metric in Eq. ͑27͒ and that in Eq. ͑28͒ is the difference between a right-handed SPC set of axes and a left-hand set.With the exception of Eqs.͑5͒-͑7͒, which are specific to the water molecule, the results in this section are general and are applicable to any rigid nonlinear molecular top.The same remains true for all the results in Sec.II B.

B. SPCs for the ellipsoid of inertia
The treatment of the three-dimensional rotations by path integral with angles is difficult for the same reasons that it is so for linear tops.It is possible to map the three-dimensional rotation of a rigid body stereographically and carry out the computation of the path integral in a "Cartesian type" set of coordinates.The simplest way is to remap the inertia ellipsoid with ⌽ 2 as follows: 3 .

͑29͒
We arrive at an expression for ⌽ 2 and its inverse by working through the quaternion parameter space associated with the three Euler angles.The detailed development of the stereographic projection map and the procedure to transform the metric tensor on the general inertia ellipsoid has been detailed elsewhere. 24,40

C. The bending coordinate
We now look for a formalism to introduce the bending DF in water.In so doing we accomplish two tasks.We are able to derive a metric tensor that allows us to perform path integral simulations of a cluster of molecules using the center of mass coordinates, the stereographic projections related to the Euler angles, and one SPC related to the bending angle.Additionally, we find a simple way to demonstrate that generally, for nonlinear tops, there are no differences in the thermodynamic averages computed from the ellipsoid of inertia space of a rigid molecule, or those computed with the infinitely stiff force constants model.
To develop the coordinate map for the space that includes the bending of water, we carry out a transformation of coordinates in two steps.First, we introduce a bending degree of freedom, which is mapped by a SPC, and we handle the rotations using Euler's angles.Second, we transform the Euler angles to SPC, ⌽ 1 maps from R 9 to the following Cartesian product of subspaces ⌽ 1 : R 9 → R 3 I 3 B 1 , where the projection coordinate 4 is the SPC of the H-O-H angle ␣ HOH : 4 spans the ring of radius r e represented by the symbol B 1 , where r e is the equilibrium length between the hydrogen atoms and the oxygen atom.We find that we can still remap I 3 with the SPCs in the usual manner. 24, , are the three Euler angles defined by the element R of the rotation group as before.It is possible to introduce a body-fixed reference configuration and rotate it, but this is a function of ␣ HOH ; the same is true for the inertia tensor: ␥ z = r e sin͑␣ HOH /2͒.͑34͒ The SPCs remapping of the trigonometric functions of ␣ HOH is obtained with double angle formulas or its additive inverse if ␣ HOH Ͼ , and is derived with geometric arguments by projecting the angular coordinate from the point on a ring of radius r e where one hydrogen is located onto a line perpendicular to the O-H vector on the opposite side of the ring ͑cf.Fig. 1͒.The expressions above are used to build the body-fixed reference frame x ͑BF͒ of Eq. ͑5͒.It can be shown that the configuration of the body-fixed frame remains in the center of mass and the inertia tensor is diagonal.With the expressions for ␣ y , ␤ y , ␥ z as a function of ␣ HOH , it is possible to obtain the mapping between the coordinates we have chosen for the subspace and the Cartesian coordinates of each of the three atoms.
We need to examine the mapping from the body-fixed to the laboratory frame.The map is expressed by Eq. ͑8͒; the elements of the first six rows of the Jacobian matrix are given by Eqs.͑13͒ and ͑14͒ for the first six rows.The last row of the Jacobian matrix of ⌽ 1 −1 is expressible as follows: R does not depend on the internal DF, only the body-fixed configuration vectors do.For the metric tensor we use the expressions in Eq. ͑16͒ for the 1 Յ Յ 3, 1Յ Յ 3 case, in Eq. ͑17͒ for the 1 Յ Յ 3, 4Յ Յ 6 case, and in Eq. ͑18͒ for the 4 Յ Յ 6, 4Յ Յ 6 case.There are three more cases to be considered: where, once more, the fact that the body-fixed frame is in the center of mass is used.The derivative operator can be moved all the way to the left since the rotation matrix does not depend on the internal coordinate.This is a general result which establishes the fact that the internal DF are also orthogonal to the center of mass coordinates.V. 4 Յ Յ 6, Ͼ 6

͑40͒
where the set ⌫ is a set of three tensors that operate in the R 3 space associated with the body-fixed Cartesian coordinates for atom i, The expressions for these are

͑44͒
To complete the metric tensor we need to explore one final case, VI.Ͼ 6, Ͼ 6, where the orthogonality of R is used.All the results produced so far in this section apply for all the internal DF.The expression in Eq. ͑45͒ does not depend on the Euler angles and r C .This result has important consequences: it implies that the submatrix part of the inverse tensor g that corresponds to the quadratic in the derivatives of the constraints as derived in Ref. 44 is independent of the Euler angles and r C .Consequently, the ratio g ͑E͒ 1/2 / g 1/2 that distinguishes the infinite spring constant limit from the curved space averages is also independent of the angles.Therefore, there are, in general, no differences between the infinite spring constant model and the inertia ellipsoid model for molecular dynamics with all the internal DF constrained.
One derives the following metric tensor for the bending water molecule with rigid stretches:

͑46͒
where the primes are used to notate the derivatives with respect to the bending coordinate.The bending coordinate is orthogonal to the rotational coordinates only because the body-fixed configuration is symmetric for every value of 4 .
In the more general case, Eq. ͑40͒ is expected to produce nonzero elements.This representation can be transformed to express the metric on the ellipsoid of inertia mapped by SPCs.It is simpler to transform directly Eq. ͑46͒ by using the map ⌽ 2 .The Jacobian matrix associated with ⌽ 2 −1 is

͑47͒
where the expressions for elements in the J block are those found in Ref. 24.It can be shown easily that J only transforms the 3 ϫ 3 rotation block of the metric tensor, and that the bending projection 4 is orthogonal to the other three SPCs:

͑48͒
where the entries of h are the usual ones, 40 and = m H ͓͑␣ y The effective mass associated with the bending projection ͑ ˙4͒ 2 is independent of the other three SPCs.g ͑E͒ 1/2 / g 1/2 does depend on the H-O-H bending DF as it is shown in Ref. 44.To see if there are any numerical differences between the infinitely stiff spring constant model and the curved space for which we derive the metric tensor g , we must find an efficient way to compute the g ͑E͒ 1/2 / g 1/2 ratio and insert the ratio into the MCPI simulation.
In Sec.II D we develop a general strategy to compute g ͑E͒ directly.

D. The bending and stretching DF
We now consider the development of a set of maps that involve all of the internal DF and the inertia ellipsoid.Mapping is performed in two stages to simplify the expression of the SPCs of the inertia ellipsoid.The maps, ⌽ 1 and ⌽ 2 , and their inverses carry out the following transformation of coordinates in R 3n : where the SPCs and the angles are the same as before and q 1 − q 3n−6 are all the internal DF ͑redundancies excluded͒.We make no distinction between the soft and the stiff internal coordinates and the particular coordinates that are used to map these.For the water molecule case specifically, q 1 is the same as 4 in the previous section, and q 2 and q 3 are the two O-H stretches r 1 and r 2 .The metric tensor is obtained by transforming Eq. ͑9͒ according to Eq. ͑10͒ in two stages, first with the Jacobian matrix of ⌽ 1 −1 and then with that of ⌽ 2 −1 .The Jacobian matrix of ⌽ 1 −1 is expressed in Eqs.͑14͒ and ͑38͒.The Jacobian matrix for ⌽ 2 −1 has a block diagonal structure, The expressions for elements in the J block are those in Refs.24 and 40.Equation ͑45͒ requires that analytical expressions for the body-fixed coordinates r i ͑BF͒ as a function of q 1 , ... ,q 3n−6 are available.
For the water molecule we need to find nine body-fixed Cartesian coordinates that are functions of 4 , r 1 , r 2 so that ͑1͒ the configuration is in the center of mass and ͑2͒ the configuration diagonalizes the moment of inertia.We begin by choosing to put the configuration in the y − z plane.Then the 1, 2 and the 1, 3 entries of the inertia tensor are zero.This choice reduces the number of degrees of freedom to six.Satisfying the first criteria produces two equations, one for the y coordinates, and one for the z coordinates; one is left with four DF.Finally, setting the 2, 3 entry of the inertia tensor to zero by rotating about the x axis produces one additional equation, reducing the number of DF to three, the same number of internal DF at our disposal.A simple way to proceed is to place the oxygen at the origin, the first hydrogen on the y axis, and the second hydrogen in the y − z plane.Then, one translates to the center of mass and rotates about the x axis until the yz moment of inertia vanishes.For the water molecule we obtain in this manner analytical expressions for the body-fixed reference configuration as a function of the internal coordinates.
It is still possible to introduce the usual SPC for the bending angle ␣ HOH , All the elements of the Jacobian matrix and the metric tensor have been represented in terms of derivatives of the rotation matrix, and the derivatives of the body-fixed coordinates.We only need to evaluate an expression for the metric tensor for values of r 1 and r 2 equal to their equilibrium value since we are only interested in fixing the distribution of 4 when we carry out simulations with the bending degree of freedom included.At r 1 = r 2 = r e we obtain a relatively simple form of the metric tensor,

͑52͒
Each of the blocks in Eq. ͑52͒ is a 3ϫ 3 block, and T denotes the transpose of .The determinant of g ͑E͒ can be easily computed with the Cholesky decomposition, g ͑E͒ = LL T , since the metric tensor is always expressible as a positive definite matrix:

E. The Riemann-Cartan curvature scalar
Both the DeWitt prepoint and the Kleinert postpoint expansion require a lattice correction, a quantum contribution of the Jacobian arising from the transformation of the Wiener measure in the manifold.The lattice correction is a scalar term, which is typically described as an effective potential, proportional to the Riemann-Cartan curvature scalar R, ͓cf.Eq. ͑2͔͒, a contraction of the Riemann curvature tensor and the ⌫ symbols are the Christoffel connections, Numerical tests show that R does not depend on the orientation, and that it is a function of the bending stereographic projection ͑or bending angle͒ only.To develop the analytic expression for R we need to choose the most convenient set of coordinates to map the R 3 I 3 B 1 manifold.Since R is scalar, it is relatively simple to use the map directly to change coordinates in our answer, without involving any elements of the Jacobian matrix.We derive the analytic expression for R using the Euler angles, as it is far simpler to find explicit equations for the connections, the inverse of the metric, and the related quantities.The first step is to find an expression for the inverse of Eq. ͑46͒ as it is needed to find the Christoffel connection coefficients in Eq. ͑56͒, and in the contraction of the Ricci tensor in Eq. ͑55͒.The block diagonal nature of g in Eq. ͑46͒ simplifies this first task, since the only nontrivial inversion takes place in the 3 ϫ 3 rotation block; only a few row elimination steps are needed to produce relatively compact expressions for the entries of g .Equation ͑56͒ yields a total of 43 nonvanishing connections for the space; however, of these only 25 expressions are distinct as the result of the symmetry property The Ricci tensor R in Eq. ͑55͒ is obtained with combinations of connection coefficients and their derivatives.In the four-dimensional subspace of interest, only 20 elements of the Riemann tensor are independent.However, even fewer are needed to compute the Ricci tensor, Additionally, the Ricci tensor is symmetric in a torsion free space, and use can be made of the antisymmetry of the last two indices, Therefore, for every element of the Ricci tensor only three elements the Riemann tensor are needed, e.g., The sum of all the Riemann tensor elements needed for the contraction gives rise to a total of 98 terms, and the simplification of the Ricci tensor elements is only marginal; a considerable number of terms for each element of R remain.Therefore, the last contraction in Eq. ͑55͒, which involves numerous accruing, factoring, and cancellations, is best left for symbolic processing software. 45

͑62͒
where all the primes denote derivatives with respect to the bending coordinate.There are several important features in this expression.͑a͒ Equation ͑62͒ does not depend on the Euler angles, neither explicitly nor implicitly.This agrees with our observation in preliminary numerical explorations.͑b͒ If we set all the derivatives with respect to the bending coordinate to zero, all but the first term vanish, and the first term is the curvature expression for the general ellipsoid of inertia of a rigid body as found in the literature.In the rigid case, the curvature is constant and has no effects on the dynamics.͑c͒ The expression in Eq. ͑62͒ is independent of the coordinate chosen to map the bending degree of freedom.

III. NUMERICAL TESTS
We have shown that it is possible, in principle, to integrate the Wiener measure on the R 3 I 3 B 1 manifold, and in a whole class of Cartesian products of these, In practice, the fact remains that these spaces are mixtures of intramolecular and intermolecular DF, which have disparate time scales.The direct use of Monte Carlo methods and the estimate of thermodynamic properties produce inefficient algorithms even when the latest numerical technology is em-ployed.The inefficiency is particularly severe in the regions of temperature where the dynamics are dominated by the motion along intermolecular DF.Nevertheless, the present development is useful to study in detail the couplings among the two sets of degrees of freedom.There are three potential sources of couplings that perturb the "separability" of intermolecular and intramolecular DF: ͑1͒ the existence of off diagonal terms in the metric tensor, ͑2͒ the existence of Christoffel connections ͑and higher de-rivatives͒ of the metric, and ͑3͒ the existence of off diagonal terms in the Hessian ͑and higher derivatives͒ of the potential.
The first two in the list are of "geometric" origin since these are specific properties of the metric tensor, whereas the last one is a property of the potential energy surface.This separation is of importance since the methods that one must employ to study these are radically different.The effects that the existence of off diagonal terms in the Hessian of the potential have on the dynamics and thermodynamics can be studied relatively well with traditional methods as, e.g., with normal mode analysis.The geometric couplings, on the other hand, are from first principles, resulting from our choice of a convenient partition of R 3n so as to "isolate" high frequency DF.
The strategy that we develop to study more closely geometric couplings is to choose a sufficiently simple system for which the existence of off diagonal terms in the Hessian ͑and higher derivatives͒ of the potential can be safely ignored.We use two separate sets of random walks in the classical and quantum limit to monitor the distribution of the low frequency DF, the first with the high frequency DF under study constrained and the second set with the high frequency DF under study sampled.The water monomer with rigid stretches, subjected to an external field, is one such example.
The potential energy model we consider has two terms, over the I 3 B 1 manifold.We omit the translation DF in the treatment since there are no couplings of any nature between the translations and the remaining DFs.We choose the interaction with the external field as a nontrivial sinusoidal function of the Euler angles, which, when expressed in terms of SPCs, takes the following form: Its shape and the value of V 0 = 0.0400 hartree are chosen to mimic approximately the interaction of a molecule of water adsorbed by a water cluster.The intramolecular interaction is simply a harmonic potential, when expressed in terms of the bending angle The parameters, b w = 0.160 988 4 hartree rad −2 and r OH,e = 1.809 23 bohr, ␣ eq = 1.823 87 radϷ 104.5°, are chosen to produce the following properties for this system.The bending frequency is 1594.3cm −1 , the experimental value in the gas phase.The harmonic approximation of the zero point energy is 0.003 63 hartree.The equilibrium value of the bending SPC is eq 4 = 2.9264 bohr.At eq 4 the inertia tensor eigenvalues are I 1 , I 2 , I 3 Ϸ 12 571, 8203, 4367 atomic units, respectively.The minimum of the overall potential is V min = −0.01732 hartree at 1 = 2 = 3 = ͱ 4 / 3, ␣ HOH = ␣ eq .Let g 1/2 be the square root of the determinant of the metric in I 3 B 1 , and g ͑E͒ 1/2 the square root of the determinant of the metric in R 9 as explained in Sec.II D. The ratio g ͑E͒ 1/2 / g 1/2 , as a function of 4 , is displayed in the graph in Fig. 2. As explained earlier, g ͑E͒ 1/2 / g 1/2 is independent of the orientation DF, a convenient feature that allows us to compute the ratio once and access its value from tabulated data during the walk.The classical Jacobian g ͑E͒ 1/2 can be computed at every point with a much cheaper evaluation of g 1/2 and a multiplication of g 1/2 by the value of g ͑E͒ 1/2 / g 1/2 interpolated linearly between its tabulated values.A similar interpolation between tabulated values is used to speed up the computation of the Riemann curvature scalar ͓cf.Eq. ͑62͔͒.With the chosen set of parameters, the curvature evaluates as shown in Fig. 3. Finally, a graph of V intra as a function of the SPC is presented in Fig. 4.
FIG. 2. The ratio g ͑E͒ 1/2 / g 1/2 as a function of the bending SPC.g ͑E͒ 1/2 is the square root of the determinant of the metric g ͑E͒ in R 9 , mapped by center of mass, orientation, and all internal DF coordinates.g 1/2 is the square root of the determinant of the metric g over the R 3 I 3 B 1 manifold.The units of the ratio are those of mass.g ͑E͒ 1/2 / g 1/2 is convenient for the fast computation of g ͑E͒ 1/2 , the classical Jacobian in simulations.
FIG. 3. The Riemann Cartan curvature scalar of the R 3 I 3 B 1 manifold as a function of the SPC 4 ͓cf.Eq. ͑62͔͒.

204107-10
To interpret the x axis of the graphs in Figs.2-4 one needs to refer back to Fig. 1.As eq 4 approaches zero, the bending angle approaches 180 deg; as eq 4 approaches +ϱ, the bending angle approaches zero.The rise in the potential energy moving to the left from equilibrium is steeper because the span from the equilibrium angle to 180 is contracted into approximately 3 bohr by the bending SPC, and is expanded greatly by the same when moving to the right of the equilibrium point.
The MCPI simulations are performed using reweighted random series ͑RRS͒.The RRS method is a finite expansion of the path with 4k m terms,

͑67͒
We use Fourier-Wiener basis, where since the RRS method 17,18 with these bases has been shown to retain its cubic convergence in non-Euclidean spaces. 24,39he expression for the importance sampling is derived from the density matrix, 24 RW ͑q,qЈ,␤͒ ϰ ͵ d͓a͔ r exp͕− ␤ ͵ 0 1 duU͑q ˜͑u͖͒͒, ͑71͒ where and where all four terms on the right-hand side are functions of q ˜͑u͒; we drop the argument of these functions to improve clarity.In this study we compute ͗V e ͘ , ͗E͘, and ͗V intra ͘. ͗V e ͘ is a good indicator of the perturbation in the sample of configuration space between simulations in I 3 B 1 and those in I 3 .
The Metropolis 46,47 algorithm is used to generate four separate walks over which values of V e are accumulated and averaged.Each walk consists of 11 million moves; the first million is used to reach the asymptotic distribution, the remainder of the walk is used to compute the averages.The entire procedure is repeated 25 times to produce a block average for ͗V e ͘ and the standard error in the mean from which to compute the error bars; these represent the 95% confidence interval.The results are presented in Figs.5-7.Incremental values of k m ͓cf.Eq. ͑67͔͒ from 0, for a classical simulations, up to 16 are used.The striking feature in Fig. 5 is the excellent agreement between the results obtained with

IV. CONCLUSIONS
In a recent article we have demonstrated that MCPI simulations of rigid molecular condensed matter, when all intramolecular DF are constrained, yield massive efficiency gains. 24However, constraining all the intramolecular DF is not always satisfactory.9][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64] In our recent work 24 on water clusters we constrain all the intramolecular DF.It is not obvious that the resulting simulation should yield satisfactory results at temperatures above 250 K, where excitation of the bending degree of freedom may have a detectable effect on sensitive thermodynamic properties.As water molecules condense, redshifting could lower this temperature by several degrees as well.The work reported in this article contains the first fundamental steps necessary to develop efficient algorithms for walks and estimators of important thermodynamic properties obtained from MCPI simulations of a cluster of molecules in which some selected internal DF are constrained by infinite spring constant, while other, softer, internal modes are included in the configuration space.We show that SPC maps can be used when the soft internal DF are bending modes.We develop a systematic method to derive the necessary expressions for the metric tensor when the space is mapped with Cartesian coordinates for the translations, stereographic projections for the orientations, and any number of non-redundant internal degrees of freedom that are not constrained.We find expressions for the correct classical and quantum Jacobians, and ways to compute these efficiently.
A careful study of the couplings that originate from the geometry of the chosen space, which is the subject of the numerical tests we report, is important for two fundamental reasons.First, it is a necessary step for the development of new theories, such as a split kinetic energy operator approach, for example, and to develop expressions for new efficient estimators.Second, the couplings among sets of DF are important when estimators for spectroscopic properties associated with high frequency DF are developed for MCPI simulations.To better grasp the relevance of the outcome of our numerical tests, let us consider the expression for the Laplace-Beltrami operator, on the I 3 B 1 manifold.The terms of the operator contain either g or its derivatives, and since the spaces R 3 I 3 B 1 are mutually orthogonal, from the block diagonal structure of the metric, there are no cross terms that mix derivatives among DF belonging to these subspaces.We symbolize this with the following expression, where the superscripts of the three terms on the right-hand side represent the space of coordinates with respect to which the derivatives are taken in the operators, while the arguments represent the functional dependence on configuration space of the operators.For example, the derivative operators for I 3 contain multiplicative factors that depend on both the orientation SPC, as well as the value of the bending coordinate.It is such dependence that gives rise to the geometric couplings we have investigated.The fact that we find no significant couplings, as evident in Fig. 5, suggests that one could drop the mixed dependence of the Laplace-Beltrami operator and to a good approximation obtain The last equation is simply a statement of the adiabatic approximation.However, it is remarkable that, in order to formulate Eq. ͑75͒, no appeal to the time scale difference ͑which is presumed to diminish the effects of dynamic couplings alone͒ is ever made.The simulations labeled "bending" in Fig. 5 are obtained by developing the propagation from the last two terms on the r.h.s. of Eq. ͑74͒.The simulations labeled "rigid" in Fig. 5 correspond to the second term on the r.h.s. of Eq. ͑75͒ only.The average external potential ͑and the difference between the two simulations͒ is a direct measure of how much the bending affects the orientation.The results in Fig. 5 can be interpreted to mean that the geometric couplings from the bending have a smaller effect than the statistical error on an average interaction felt by a molecule of water on the surface of a cluster.Obviously, at some level such coupling can be measured, and a good indicator would be the values of the Christoffel coefficients that connect bending and orientations.The connecting Christoffel coefficients are numerous, and not all are vanishing, but they are clearly small.Similar arguments apply for dynamic couplings as well.Off diagonal elements of Hessian matrices are never zero, but they might be smaller than the statistical error on the average potential energy coming from the relatively larger diagonal elements.In Euclidean spaces, that evidence is sufficient to justify splitting the Hamiltonian and regard the adiabatic approximation as excellent.In non-Euclidean ͑or even curvilinear͒ spaces, the present article argues that other sources of coupling exist, that these couplings arise from the properties of the metric tensor, and that they need to be researched carefully.The bending mode is, from the path integral prospective, a "new" type of space which must not be mistaken with the typical mono-dimensional rotation.The element of the metric tensor is constructed from the masses of the three atoms involved in the bending, and the derivatives of the bodyfixed frame coordinates of these.The dependence on the space is implicit and, therefore, the expression for g 77 in Eq. ͑46͒ is independent of the coordinate chosen to map B 1 .However, the most striking difference between the bending DF and a typical mono-dimensional rotation can be observed by considering the expression for ⌬ LB where g I 3 is the determinant of the metric tensor block associated with the rotation space, and the primes notate derivatives with respect to the bending DF.This expression is also independent of the coordinate chosen to map B 1 , remains nearly unchanged in the adiabatic limit, and is quite different from the Laplace-Beltrami operator for the particle in a ring.These are the reasons that compel us to call the bending space B 1 , and not S 1 , which should be reserved for the monodimensional rotation.
The next important step necessary to achieve our objectives is a careful study, perhaps from first principle, of the dynamic couplings in water clusters.Unlike the geometric couplings computed in the present paper, the dynamic couplings require a carefully constructed potential energy surface.Many flexible extensions of the rigid water PES models ͑e.g., the TIP4P or the SCPE͒ exist.However, the dynamic couplings we seek require additional theoretical or experimental studies.
The results obtained in the present investigation, namely the orthogonality among the rotation and the bending, and the lack of geometric couplings can be generalized only to triatomic molecules with C 2v symmetry.As stated in Sec.II C, the bending coordinate is orthogonal to the rotational coordinates only because the body-fixed configuration is symmetric for every value of 4 .Bending DF may not be orthogonal to rotations in less symmetric molecules.Nonorthogonality gives rise to distinct types of coupling terms whose importance must be carefully scrutinized with the theoretical and numerical methods we employ here.
The method we develop here can be clearly generalized to study the differential geometry with any kind of relatively low frequency internal degree of freedom.Equation ͑45͒, for example, is quite general, as it applies to any kind of internal degree of freedom.However, we caution the reader that, while the procedure yields the correct expressions for the metric tensors and classical and quantum Jacobians, the fact that these exist and can be obtained does not assure that the resulting path integral converges.The Wiener measure in manifolds remains a rich subject, and the theoretical methods that the community have developed are far from being black boxes.As a clear example of how partitioning a space, along with certain choices of coordinates, can lead to nonconvergent path integrals, one need only consider the well known problem of path collapsing along the radial degree of a rotating linear top: When mapping the rotations with traditional polar angles, a purely geometric term in the radial path integral for the l = 0 state is infinitely attractive as r approaches zero. 29Special remappings of the radial manifold have been employed to handle the path collapse problem for a rotating linear top. 29On the other hand, the procedure we develop here can be used to study more general types of bending, and torsional manifolds.It is likely that all these can be mapped with SPCs, and that these coordinates are orthogonal to the rotation DF.

ACKNOWLEDGMENTS
This work has been supported by the National Science Foundation ͑Grant No. CHE0554922͒.Additionally, E.C. acknowledges the donors of the Petroleum Research Fund, administered by the ACS ͑Grant No. 40946-B6͒, The Stacy Ann Vitetta '82 Professorship Fund, and The Ellington Beavers Fund for Intellectual Inquiry from Arcadia University for partial support of this research.E.C. would also like to thank fellow faculty members in the Mathematics Department at Arcadia University, for access to symbolic processing software.

APPENDIX: METRIC TENSOR FROM PASSIVE ROTATION OF A RIGID WATER MOLECULE
It is straightforward to show that for the passive rotation case the relevant R 3 tensors, obtained by taking derivatives of the transpose of the matrix in Eq. ͑4͒, are 2 + cos 2 cos sin ͑1 − cos 2 ͒ − sin cos sin cos sin ͑1 − cos 2 ͒ cos 2 cos 2 + sin 2 cos sin cos − sin cos sin cos sin cos sin 2 .

FIG. 1 .
FIG.1.A sketch of the geometric definition of the two coordinates used to map the bending manifold B 1 .The stereographic projection is obtained by finding the intersection point between the line that originates at H, passes through the point on the circle at ␣ HOH , and the tangent line that lays perpendicular to the line between the center and H.

FIG. 4 .FIG. 6 .
FIG. 4. Graph of V intra as a function of the SPC 4 .FIG. 5. Value of ͗V e ͘ in the I 3 manifold ͓classical ͑white circles͒, quantum ͑white squares͒ labeled as "rigid"͔ compared to the values of ͗V e ͘ computed with the I 3 B 1 manifold ͓classical ‫,͒ء͑‬ quantum ͑+͔͒.The quantum values are obtained with k m = 4 in both cases.
and bending J. Chem.Phys.128, 204107 ͑2008͒ ⌫ 55 Ј = cos 2 sin 2 + cos 2 cos sin ͑cos 2 − 1͒ − sin cos sin cos sin ͑cos 2 − 1͒ cos 2 cos 2 + sin 2 − cos sin cos − sin cos sin − m t , g 44 = I 11 cos 2 + I 22 sin 2 , g 46 = ͑I 11 − I 22 ͒sin cos sin , 2 66 = I 11 sin 2 sin 2 + I 22 sin 2 cos 2 + I 33 cos2.͑27͒Equation͑27͒ is one of the few equivalent forms for the metric tensor in the inertia ellipsoid mapped by Euler angles; the expression applies in general to n particles treated as a single rigid body.We can develop an equivalent form of the metric if we handle the rotation passively ͑i.e., the rotation is applied to the body-fixed axis, rather than the coordinates of the atoms͒.The transformation map in the passive case is Finally, since some elements of the inverse of the metric tensor are zero, a number of elements of the Ricci tensor can be omitted.Only R 44 , R 45 , R 46 , R 55 , R 56 , R 66 , and R 77 are needed, for a total of 21 elements of the Riemann tensor.