Full Dimensional (15D) QuantumDynamical Simulation of the Protonated WaterDimer I: Hamiltonian Setup and Analysis of the Ground Vibrational State
Abstract
Quantumdynamical fulldimensional (15D) calculations are reported for the protonated water dimer (HO) using the multiconfiguration timedependent Hartree (MCTDH) method. The dynamics is described by curvilinear coordinates. The expression of the kinetic energy operator in this set of coordinates is given and its derivation, following the polyspherical method, is discussed. The PES employed is that of Huang et al. [JCP, 122, 044308, (2005)]. A scheme for the representation of the potential energy surface (PES) is discussed which is based on a high dimensional model representation scheme (cutHDMR), but modified to take advantage of the modecombination representation of the vibrational wavefunction used in MCTDH. The convergence of the PES expansion used is quantified and evidence is provided that it correctly reproduces the reference PES at least for the range of energies of interest. The reported zero point energy of the system is converged with respect to the MCTDH expansion and in excellent agreement ( cm below) with the diffusion Monte Carlo result on the PES of Huang et al. The highly fluxional nature of the cation is accounted for through use of curvilinear coordinates. The system is found to interconvert between equivalent minima through wagging and internal rotation motions already when in the ground vibrationalstate, i.e., T=0. It is shown that a converged quantumdynamical description of such a flexible, multiminima system is possible.
I Introduction
The protonated waterdimer HO (Zundel cation) is the smallest system in which an excess proton is shared between two water molecules. Much effort has been devoted to this system due to the importance that the solvated proton has in several areas of chemistry and biology. In recent years a fast development of the spectroscopical techniques available to probe the vibrational dynamics of ionic species in the gas phase has taken place, and several studies have appeared around the HO system Asmis et al. (2003); Fridgen et al. (2004); Headrick et al. (2004); Hammer et al. (2005) and other more complex clusters and molecules Headrick et al. (2005); Roscioli et al. (2007). In order to achieve a satisfactory understanding of the spectroscopy and dynamics accurate theoretical results are needed in parallel and several works have appeared to fill this gap which are based on different theoretical approaches, from classical to quantumdynamical methods Vener et al. (2001); Dai et al. (2003); Huang et al. (2005); Hammer et al. (2005); Sauer and Dobler (2005); Kaledin et al. (2006). Accurate quantumdynamical simulations of the dynamics and spectroscopy of the system require of an accurate reference potential energy surface (PES). Very accurate PES and dipole moment surfaces (DMS) have been recently produced by Huang et al. Huang et al. (2005) for HO. They are based on several tenthousands of coupledcluster calculations combined with a clever fitting algorithm which is based on a redundant set of coordinates, namely all atomatom distances. The PES is able to describe the floppy motions occurring at typical energies of excitation in the linear IR regime and dissociates correctly Huang et al. (2005). Several works have already appeared in which the PES and DMS of Huang et al. Huang et al. (2005) were used McCoy et al. (2005); Hammer et al. (2005); Kaledin et al. (2006). In Ref. McCoy et al. (2005) fulldimensional vibrational calculations for HO were undertaken using diffusion Monte Carlo (DMC) and variational wavefunction methods. The zero point energy (ZPE), some vibrational excitedstate energies and properties of the ground vibrational state of the system are reported and discussed. The PES of Huang et al. was also used in Ref. Hammer et al. (2005) where the experimental IR spectrum of the system was reported together with the calculation of relevant excitedvibrational states.
As will be later discussed, HO features several symmetryequivalent minimum energy structures and large amplitude motions between different regions connected by low energy barriers. For such systems curvilinear coordinates, e.g. bondbond angle, dihedral angle, or internal rotation coordinates, offer an advantage over rectilinear coordinates since they provide a physically meaningful description of the different largeamplitude molecular motions. As a remark, an attempt was made to describe the system by a set of rectilinear coordinates, since they result in a simple expression of the kinetic energy operator (KEO). The resulting Hamiltonian was unsuccessful at describing of variuos largeamplitude displacements, which resulted in abandonement of such an approach and introduction of a curvilinearcoordinates based Hamiltonian. When deriving the KEO for a Hamitonian in curvilinear coordinates one may define those in terms of rectilinear atom motions and use standard differential calculus to obtain the desired expression Podolsky (1928). This approach becomes extremely complex and error prone already for systems of a few atoms. Another possibility is to use an approximate KEO where some crossterms are neglected, thus simplifying the derivation. In this case, however, an uncontrolled source of error is introduced in the calculation. These drawbacks are overcome by employing the polyspherical approach Gatti et al. (1998a); Gatti (1999); Gatti et al. (2001) when defining the coordinate set and deriving the corresponding KEO. In the polyspherical method the system is described by a set of vectors of any kind (e.g. Jacobi, valence, …). The kinetic energy is given in terms of the variation in length and orientation of these vectors, the latter is defined by angles with respect to a body fixed frame. Following the polyspherical method the KEO of the system at hand is derived in a systematic way without use of differential calculus. The polyspherical approach has already been discussed in a general framework as well as for some molecular systems (see for instance Gatti et al. (1998a); Gatti (1999); Gatti et al. (2001); Gatti and Iung (2003)). In this paper the full derivation of the KEO in a set of polyspherical coordinates is illustrated for the HO cation and all the key steps are carefully discussed. The derivation is followed stepbystep for this rather complex system, thus providing the basic guidelines to be followed for the derivation of KEOs for complex molecular systems and clusters.
In quantumdynamical wavefunctionbased studies, as the one being introduced here, it is convenient to represent the operators on a discrete variable representation (DVR) grid, which in the multidimensional case is the direct product of onedimensional DVRs. The potential operator is then given by the value of the PES at each point of the grid. In the present case, however, the resulting primitive grid has more than points. This number makes clear that the potential must be represented in a more compact form to make the calculations feasible. Algorithms exist which take advantage of the fact that correlation usually involves the concerted evolution of only a small number of coordinates as compared to the total number. Hierarchical representations of the PES are then constructed including potentialfunction terms up to a given number of coordinates Bowman et al. (2003); Rabitz and Alis (1999); Li et al. (2001). We discuss here a variation of a hierarchical representation of the PES in which the coordinates are grouped together into modes, and the modes are the basic units that define the hierarchical expansion. It will be shown that this approach allows for a fast convergence of the PES representation if the modes are defined as groups of the most strongly correlated coordinates.
The present work provides new reference results on the properties of the ground vibrational state of HO. The convergence of the given results is established by comparison to the DMC results McCoy et al. (2005) on the reference PES of Huang et al. Huang et al. (2005). These results, together with the Hamiltonian defined here, set also the theoretical and methodological framework used in the companion paper Vendrell et al. (2007) where the IR spectrum and vibrational dynamics of HO are analyzed.
The paper is organized as follows. Section II presents a brief description of the multiconfiguration timedependent Hartree (MCTDH) method. Section III details the construction of the Hamiltonian operator for the HO system. The derivation of the KEO is discussed in III.1, while the construction of the potential is detailed in III.2. Section IV.1 discusses the calculation of the ground vibrational state of the system and gives a comparison to available results. In section IV.2 some properties of the system in relation to its ground vibrationalstate wavefunction are discussed. In section IV.3 the quality of the potential expansion is analyzed and discussed. Section V provides some general conclusions.
Ii Brief description of the MCTDH method
The Multiconfiguration timedependent Hartree (MCTDH) method Meyer et al. (1990); Manthe et al. (1992); Beck et al. (2000); Meyer and Worth (2003) is a general algorithm to solve the timedependent Schrödinger equation. The MCTDH wave function is expanded in a sum of products of so–called single–particle functions (SPFs). The SPFs, , may be one– or multi–dimensional functions and, in the latter case, the coordinate is a collective one, . As the SPFs are time–dependent, they follow the wave packet and often a rather small number of SPFs suffices for convergence.
The ansatz for the MCTDH wavefunction reads
(1)  
where denotes the number of degrees of freedom and the number of MCTDH particles, also called combined modes. There are SPFs for the ’th particle. The denote the MCTDH expansion coefficients and the configurations, or Hartreeproducts, are products of SPFs, implicitly defined by Eq. (1). The SPFs are finally represented by linear combinations of timeindependent primitive basis functions or DVR grids.
The MCTDH equations of motion are derived by applying the DiracFrenkel variational principle to the ansatz Eq. (1). After some algebra, one obtains
(2)  
(3) 
where a vector notation, , is used. Details on the derivation, the definitions of the meanfield , the density and the projector , as well as more general results, can be found in Refs. Manthe et al. (1992); Beck et al. (2000); Meyer and Worth (2003). Here and in the following (except for Section III) we use a unit system with .
The MCTDH equations conserve the norm and, for timeindependent Hamiltonians, the total energy. MCTDH simplifies to TimeDependent Hartree when setting all . Increasing the recovers more and more correlation, until finally, when equals the number of primitive functions, the standard method (i. e. propagating the wave packet on the primitive basis) is used. It is important to note, that MCTDH uses variationally optimal SPFs, because this ensures early convergence.
The solution of the equations of motion requires to build the mean–fields at every time–step. A fast evaluation of the mean–fields is hence essential. To this end the Hamiltonian is written as a sum of products of mono–particle operators:
(4) 
where operates on the th particle only and where the are numbers. In this case the matrixelements of the Hamiltonian can be expressed by a sum of products of mono–mode integrals, the evaluation of which is fast.
(5) 
Similar expressions apply to the matrix of meanfields . For further information on the MCTDH method, see http://www.pci.uniheidelberg.de/tc/usr/mctdh/.
Iii System Hamiltonian
iii.1 Kinetic energy operator for HO
The derivation of the exact
kinetic energy operator
for HO is based
on a polyspherical approach
which has been devised in previous articles
Gatti et al. (1998a, b); Iung et al. (1999); Gatti (1999); Gatti et al. (2001); Gatti and Iung (2003).
This approach can be seen as a
very efficient way to obtain
kinetic energy operators for the
family of polyspherical coordinates.
The formalism can be applied
whatever the number of atoms
and whatever the set of vectors: Jacobi,
Radau, valence, satellite, etc.
The formalism is not restricted to total J=0
and hence the operator may include overall rotation
and Coriolis coupling.
The protonated water dimer system is described by six Jacobi vectors. The choice of the Jacobi vectors is not unique and several clustering schemes are possible, one natural choice for HO is given in Figure 1.
FIGURE 1 AROUND HERE
The polyspherical approach straightforwardly provides the expression of the kinetic energy operator in terms of the angular momenta associated with the vectors describing the system. Here we use the technique of “separation into two subsystems” which is described in Ref. Gatti, 1999 (see Eq. (37) there). The full kinetic energy operator reads (we use the notation and throughout the paper):
with

R, R, R, R, R, r, are the lengths of the vectors , , , , , , respectively.

= , = , = , and .

is the angular momentum associated with .

is the angular momentum associated with and (= ) is the total angular momentum associated with monomer A (B).

is the angular momentum associated with the proton.

is the total angular momentum of the system. = + + + and is the angular momentum associated with .

E2 is the frame resulting from the two first Euler rotations starting from the space fixed frame, such that the z axis lies parallel to . In other words, the two first Euler angles and are identical to the two spherical angles of in the space fixed frame.

BFA (B) is the body fixed frame associated with monomer A (B). The two first Euler angles and of monomer A (B) are defined such as z lies parallel to (in other words they are identical to the two spherical angles of in the E2 frame). The third Euler angle is defined such as remains parallel to the (, ) half plane. This definition of BFA(B) is similar to the one recently used to describe the water dimer (see e.g. Ref Leforestier et al. (2002)). Note however that we have changed the definition of the z axis which is now parallel to the vector between the two hydrogen atoms and not to the vector joining the center of mass of H to the oxygen atom as in Ref. Leforestier et al. (2002). This change was performed to avoid the singularity which appears in the KEO when and z are parallel.
It should be emphasized that all the angular momenta appearing
in the kinetic energy operator are all computed in the space
fixed frame but projected onto the axes of different frames.
This last point, as well as the properties of
the corresponding projections,
is addressed in detail in Ref. Gatti and Nauts (2003).
We recall here only the main point, i.e. the fact that if a
vector is not involved in the definition of a frame F, the
expression of the projections of the corresponding
angular momentum onto the Faxes in terms of
the coordinates in this frame, is identical to the usual one in
a space fixed frame.
This very useful property will be
utilized several times in the following
for instance to obtain Eqs.
(13,14,15,25,28,29).
We now slightly change the volume element for the lengths of the vectors (as usual), and assume J=0. The operator can then be recast in the following form:
which is to be used with the volume element :
(8) 
with
We have 16 degrees of freedom instead of 15 since we have not defined the third Euler angle yet. The angles are depicted in Figure 2:
FIGURE 2 AROUND HERE

and are the spherical angles of the proton in the E2 frame: is the angle between z and the vector , describes the rotation of around z.

and are the spherical angles of in the E2 frame. The angle describes the rocking motion of water A fragment while describes the rotation of the water A fragment around the axis (the same for B).

and are the spherical angles of in the E2A frame, the E2A frame is the frame obtained after the two Euler rotations and with respect to the E2 frame. Consequently, is the angle between and and describes the rotation of around . Hence describes the wagging motion of the water A fragment (the same for B). The fact that we have used the separation into two subsystems for the two water monomers (see Eq. (37) in Ref. Gatti (1999)) explains why the angles of are defined with respect to and not with respect to as in the standard formulation of the polyspherical approach (see Eq. (65) in Ref. Gatti et al. (1998a)). This allows us to have purely intramonomer angles instead of angles mixing the intramonomer and intermonomer motions.
Since = (here we can use rather than since the projections of onto the axes of the E2 or BFA (B) frames are hermitian: see Eqs. (13,15,25,29) below), we can rewrite the operator as:
This operator could be straightforwardly used along with an adequate primitive basis set of spherical harmonics and Wigner rotationmatrix elements. However, in this work, we prefer to employ a direct product primitive basis such as a DVR for all the degrees of freedom (only because the latter basis set is numerically more efficient for our problem). We have then to explicitly express of the angular momenta in terms of the coordinates, i.e. we have to detail the following terms:

(1) , , ,

(2) ,

(3)

(4)

(5)
but before doing so we have to further specify the coordinates:

(i) to avoid the singularity we use the BF Cartesian coordinates , , and for the proton (the definition of BF frame is given in (ii)). Hence we use, .

(ii) we slightly change the coordinate system (from E2 to BF) and explicitly define the third Euler angle of the system and then the BF frame. We define such that it is identical to the first Euler angle of monomer A (here, we slightly break the symmetry between the two monomers). Consequently, we have:
since J=0, we can then put:
(11) and find
Note that this transformation does not affect the other angles.

(iii) in order to have hermitian conjugate momenta for the angles, we use u = , u = , u = , u = .
For each angle , we have (with u = ): and .
We have now 15 degrees of freedom:
Let us now detail all the terms appearing in the kinetic energy operator:

(1) For , the situation is simple since monomer B is not involved in the definition of the body fixed frame of the whole system. Thus, we have :
(13) and
(14) Similarly can be seen as the total angular momentum of monomer A projected onto the axes of the Body fixed frame of the monomer. We obtain the usual expression:
(15) and
(16) However, since monomer A is involved in the definition of the third Euler angle of the system, we apply the change of coordinates, Eqs. (III.1) and (III.1):
For the projections of onto to the BFA (B) axes, the situation is not straightforward since is involved in the definition of BFA (B). However, since we have used the same convention for the BF frames as in our previous articles, the expressions of these projections is already known (identical to Eq. (47a) in Ref. Gatti et al. (1998a) for instance):
(18) (This expression can be obtained starting from the E2A frame which is the frame resulting from the two first Euler rotations and . The expression of the projections of onto this E2A frame are the usual ones since is not involved in the definition of this frame. Applying then the third Euler rotation and the coordinate transformation leading to the BFA coordinates yields Eq. (18)). Note that is not hermitian (this is due to the fact that is involved in the definition of BFA). However, at the end we obtain the usual expression:
(19) and exactly the same for B:
(20) 
(2) The terms in are not obvious since is involved in the definition of BFA (B). First, let us explicitly highlight the hermiticity of the term:

(3) For the term in , we note that the expression of the projections of onto the E2 can be seen as those of total angular momentum of a system onto the axes of the Space fixed frame. Indeed, since monomers A and B are not involved in the definition of E2, we have (similar to Eq. (18) in Ref. Gatti et al. (1998a)):