A diabatic parameterization of the twofold ground state potential energy surface of the H2O-OH molecular complex

A diabatic parameterization of the twofold ground state potential energy surface of the H2O-OH molecular complex E. Galbis, E. Giglio, and B. Gervais Citation: The Journal of Chemical Physics 139, 164313 (2013); doi: 10.1063/1.4826319 View online: http://dx.doi.org/10.1063/1.4826319 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/16?ver=pdfcov Published by the AIP Publishing Articles you may be interested in A global ab initio potential energy surface for the X 2 A ′ ground state of the Si + OH → SiO + H reaction J. Chem. Phys. 139, 204305 (2013); 10.1063/1.4832324 Correlated ab initio investigations on the intermolecular and intramolecular potential energy surfaces in the ground electronic state of the O 2 − ( X Π g 2 ) − HF ( X Σ + 1 ) complex J. Chem. Phys. 138, 014304 (2013); 10.1063/1.4772653 Ab initio and analytical intermolecular potential for Cl O – H 2 O J. Chem. Phys. 126, 114304 (2007); 10.1063/1.2566537 An ab initio correlated study of the potential energy surface for the HOBr.H 2 O complex J. Chem. Phys. 121, 141 (2004); 10.1063/1.1755191 Potential energy surface and rovibrational states of the ground Ar–HI complex J. Chem. Phys. 120, 6471 (2004); 10.1063/1.1665467 Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44 THE JOURNAL OF CHEMICAL PHYSICS 139, 164313 (2013) A diabatic parameterization of the twofold ground state potential energy surface of the H2O-OH molecular complex E. Galbis, E. Giglio, and B. Gervaisa) CIMAP, unité mixte CEA-CNRS-ENSICAEN-UCBN 6252 BP 5133, F-14070 Caen, Cedex 05, France (Received 24 June 2013; accepted 7 October 2013; published online 28 October 2013) We present a matrix functional form to fit the nearly degenerated potential energy surface of the H2O-OH molecular complex. The functional form is based on second order perturbation theory, which allows us to define two diabatic states coupled together in the field of the surrounding water molecules. The fit reproduces faithfully the fine details of the potential energy surface (PES) like the crossings and the shallow barrier between the main and secondary minima. The explicit dependence of the model on polarization ensures its transferability to systems made of several water molecules. The potential is used to investigate the structural properties of the OH radical in solution by Monte Carlo simulation. The twin surface fit shows that the second PES is shifted above the ground state by typically 1600 cm−1for the configurations explored at a temperature of 300 K and a density of 1.0 g/cm3. The second PES has thus little influence on the structuring of water around the OH radical at such a temperature and density. Our study confirms that under these thermodynamic conditions, OH is a weak hydrogen acceptor. © 2013 AIP Publishing LLC.[http://dx.doi.org/10.1063/1.4826319] I. INTRODUCTION The OH radical is one of the most abundant chemical species produced in water radiolysis, with hydronium and solvated electron.1,2It is a highly reactive species, leading to oxidation reactions as a basic mechanism observed in atmospheric reaction, in radiobiology and radiotherapy, or in corrosion, for example. To investigate the underlying mechanisms, molecular dynamics is a powerful tool, which provides detailed structural and dynamical information. The necessary input for such an investigation is of course an accurate model of the potential energy surface (PES) of the system. The Carr-Parinello Molecular Dynamics (CPMD) approach uses a global PES obtained by intensive density functional theory (DFT) calculations. It has been used recently to investigate the diffusion mechanism of OH through liquid water.3Such an approach is however computationally cumbersome and thus limited regarding the sample size, propagation time, and number of trajectories. It has been shown that self-interaction correction (SIC) to DFT is mandatory to produce meaningful results, making the computation even more demanding.4 To overcome these limitations, a parameterized model PES could be used. Such an approach has been widely used to investigate material properties, and in particular those of water, by means of conventional MD or by Monte Carlo (MC) sampling of the phase space.5Such a method was used to investigate the solvation of OH in bulk and at the surface of water samples, but it was based on a rudimentary PES model.6 The H2O-OH interaction is dominated by electrostatics in the long range domain and by hydrogen bond formation around equilibrium distances. Previous works have been performed to identify the equilibrium geometries and harmonic a)Author to whom correspondence should be addressed. Electronic mail: gerv[email protected] frequencies of the complex for the main and the secondary minimum.7,8Several recent calculations performed with different computational methods confirm the energy and geometry of these minima associated with hydrogen bond formation as donor or acceptor.9,10 Experiments were also performed to analyze finely the shape of the coupled PES around the main minimum.11 The transition state associated with hydrogen atom transfer from H2O to OH was also identified.10,12–14 All these studies show that the transition gap is located relatively far above the two minima. The rest of the surface, i.e., most of the five-dimensional space associated with the geometry of the complex for rigid H2O and OH remains largely unexplored. A previous attempt to obtain an analytical representation of the PES9is based on an adiabatic representation of the PES and it was fitted on a limited set of data. The coupling between the two PES is simply ignored in such an approach. However, for an accurate representation of the PES, these couplings have to be taken into account, as we shall see. The existing fit is definitely accurate near the minima and can be used to investigate the stationary points of the PES, but this might not be sufficient to investigate the dynamics of the system, when its internal energy increases, as it is the case for condensed phase at high temperature or pressure and also in gas phase. To represent accurately the PES of the H2O-OH complex, the specificity of the OH radical should be taken into account. Besides the Kramers degeneracy, the X2ground state of OH radical is doubly degenerate. The two states differ by the occupation of the πorbitals, and thus by the orientation of the quadrupole. This feature does not exist in a restricted mean field approach and a minimal approach based on two electronic configurations is necessary to restore this difference. The aim of this work is to develop a matrix formulation of the potential energy, which is able to take into account accurately the complexity of the coupled ground states PES of 0021-9606/2013/139(16)/164313/10/$30.00 © 2013 AIP Publishing LLC139, 164313-1 Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44 164313-2 Galbis, Giglio, and Gervais J. Chem. Phys. 139, 164313 (2013) the H2O-OH complex. Such kind of analytic PES has a broad range of applications. Besides standard molecular dynamics studies of systems made of hydroxyl embedded in water, it offers the possibility to investigate specifically the role of the PES coupling by spin-orbit or non-adiabatic interactions, and it allows to study the interaction of hydroxyl with other radicals in aqueous solution. In connection with H4O+ 2PES, it will afford the investigation of relaxation dynamics of water clusters following ionization. We use the developed potential to investigate the structure of liquid water solution embedding an OH radical by means of MC sampling of the configuration space. The radial distribution function around the hydroxyl is probably difficult to obtain experimentally because of the high reaction rate between two hydroxyl radicals. To our best knowledge, only a few simulations have been performed to obtain it. Previous works in this field were based on very simple H2O-OH potentials6on the one hand, and on much more time consuming self-interaction-corrected density functional method,4on the other hand. These very different approaches lead to a quite different picture of the OH radical solvation in liquid water. Our approach provides a third alternative, realizing a trade-off between computational efficiency and accuracy. Atomic units (¯=e=me=1) are used throughout the paper, unless otherwise stated. II. REFERENCE POINTS GENERATION We performed multireference configuration interaction (MRCI) calculations with Davidson correction in order to obtain the PES for the two lowest eigenstates of the H2O-OH molecular complex. We used two reference wave functions generated from a complete active space self-consistent field (CASSCF) built on 15 electrons in 11 orbitals. To obtain a balanced description of the two reference states, we took an averaged mean field over the two lowest CASSCF states. All calculations were performed using the MOLCAS quantum chemistry code.15 We employed a small natural atomic orbital (ANO-S) basis set included in MOLCAS code for both O and H. The contraction is (10s6p3d)/[7s6p3d] for oxygen (7s3p)/[4s3p] and for hydrogen.16 The basis set superposition error (BSSE) correction was not taken into account. Our calculations of equilibrium distance Reand depth De for the main minimum are reported in Table I, where they are compared to some values collected in the literature. Our MRCI (ANO-S) compares quite well with previous MRCI calculation.8The geometries of both calculations are actually slightly different and the ANO-S basis set is slightly larger than the basis set used by Xie and Schaefer.8Our total energies are thus slightly deeper, though Deis slightly smaller than the result of Xie and Schaefer. For comparison, we computed Reand Dewith the coupled cluster with single and double excitation and peturbative triple excitation correction (CCSD(T)), which is more correlated and size-consistent. Using the CCSD(T) method with (ANO-S) basis set gives a slightly smoother potential well than our MRCI (ANO-S) calculation. Figure 1shows a scan of the PES for OH acting as a donor at several levels of approximation. The scan was per- TABLE I. Reand Defor the different approach used to scan the PES for the H2O-OH system. The values of Ref. 8are given for a bending angle θA=30◦. ReDeDe Method (a.u.) (10−3a.u.) (kcal/mol) MRCI (ANO-S) 5.50 8.92 5.59 CCSD(T) (ANO-S) 5.54 8.50 5.33 CCSD(T) (aug-cc-pVQZ) 5.46 9.15 5.74 CCSD(T) (aug-cc-pVQZ) BSSE cor. 5.48 8.82 5.53 CCSD(T) (aug-cc-pVTZ) Ref. 95.48 9.33 5.85 MRCI (TZ2P) Ref. 85.51 9.07 5.69 formed at a bending angle θA=36◦, which corresponds to the geometry of the main minimum obtained by Du et al.9We see that our MRCI (ANO-S) approach compares fairly well with the CCSD(T) calculation performed with the much larger aug-cc-VQZP basis set and taking into account the BSSE correction by means of the counterpoise method17 for each point. We are thus confident in the quality of the PES generated by using the MRCI (ANO-S) method. The H2O-OH system has nine internal degrees of freedom, corresponding to two interatomic distances for the water molecule, one interatomic distance for the hydroxyl, one bending angle for the H2O, one intermolecular distance between H2O and OH (ROO) and four angular coordinates (θA, φA,θBand φB) defining the relative orientation of the two molecules as illustrated in Figure 1(seealsoRef.18). We fixed the interatomic distances at a common values for the water molecule and the OH radical, rOH =1.84 a.u. and αHOH =104.0◦.ForH 2O-OH complex, the intramolecular degrees of freedom are well decoupled from the intermolecular degrees of freedom.8Moreover, for OH embedded in water at thermal energies, the internal vibrational modes have limited amplitudes. It is thus reasonable to fit separately the intermolecular potential energy. In order to scan correctly the PES, we vary ROO,θA,φA,θB, and φB(see Ref. 18) to obtain approximately 20 000 points for each PES. We employed 20 ROO values in the range [4.0, 12.5] a.u. with a larger density of points at short distance, where the potential energy varies quickly with ROO. III. FIT FUNCTIONAL FORM We choose conventionally x2>y2for the 2xstate and the reverse for the 2ystate. The intrinsic degeneracy between these states is lifted when an external electric field is applied on the molecule, as it is the case when OH forms a complex with additional molecules. Following Spelsberg,19 we start our analysis from perturbation theory and seek for a matrix representation of the coupled PES for the two lowest energy states of the H2O-OH molecular complex. Denoting H2O by A and OH by B, the potential energy matrix reads U(A, B)=Uxx Uxy Uyx Uyy .(1) Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44 164313-3 Galbis, Giglio, and Gervais J. Chem. Phys. 139, 164313 (2013) 5678 9 10 11 12 ROO (a.u.) -0.01 -0.008 -0.006 -0.004 -0.002 0 E (a.u.) MRCI CCSD(T) CCSD(T) (BSSE corrected) Fit z FIG. 1. Potential energy curves of H2O-OH system for different levels of calculation. For a system made of several H2O molecules, we distinguish additive and non-additive contributions. Here, the additivity refers to matrix addition, i.e., the matrix elements are added together to form the resulting many-body interaction matrix U. The additive contributions are made of dispersion, electrostatic, and short range interactions as detailed below. Among them, only the dispersion was not given an off-diagonal term and is thus additive in the usual sense of pair energy additivity. The non-additive contribution comes from induction only. For additive contributions we thus evaluate the matrix elements for each pair (A,B)toformV(B) =AV(A, B). For induction contribution, the field created by the surrounding molecule Con each molecule Dis first evaluated, and then the matrix elements evaluated to form the induction energy matrix W, as we shall see below. The energy of one OH radical interacting with several water molecules is thus not simply the sum of the pair interaction energy of OH with each water molecule. It results from the diagonalization of the matrix U=V+W. For a system made of two OH radicals, the size extensivity would be obtained by increasing the rank of the matrix to take into account the 4 states asymptotically degenerated for each spin state. In such a case, a model PES could be built, like in the empirical valence bond approach,20 by coupling the two pairs of states. The coupling needs however to be determined from additional quantum chemistry computations. To evaluate the interaction matrix of OH with the surrounding water molecules, we first choose a reference frame to orient the OH molecule and to define xand y. The matrix U is real and symmetric, so that the energies of the system, computed as the eigenvalues of U, are of course real. The energies of the system should also be independent on the frame chosen to orient the OH states 2xand 2yand invariant under the exchange of xand y.Ifφis the azimuthal angle defining the orientation of 2xin the reference frame, the symmetry constraints can be expressed as Uxx(φ)=U0+U2cos(2φ),(2) Uyy(φ)=U0−U2cos(2φ),(3) Uxy(φ)=U2sin(2φ).(4) Both U0and U2depend on the coordinates that express the relative positions of the molecules. U0is the part invariant under rotation around the OH molecular axis. It contains, for example, the dipole type matrix element common to both 2x and 2ystates. U2is related to the difference between these two diabatic states and contains, for example, the asymptotic quadrupole contribution from OH, which is different for the 2xand 2ystates. From an analysis based on second order perturbation theory, the long range interaction energy can be decomposed into electrostatic, polarization, and dispersion contributions. Each contribution can be determined from individual properties of each molecules.19 In the short range limit, such an analysis is no longer applicable and hydrogen bonding becomes important. We thus damp the electrostatic, polarization, and dispersion interactions and introduce a short range contribution. The equivalence between the two diabatic states is however maintained, and the functional form of Uxx and Uyy and Uxy uses the above decomposition. Both molecules carry static multipoles and the electrostatic contribution dominates by far the long range interaction. Like many authors in this field,21–23 we give to each atom a set of multipoles. The matrix elements are expressed as a function Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44 164313-4 Galbis, Giglio, and Gervais J. Chem. Phys. 139, 164313 (2013) of the atomic multipoles as follows: Vxy = a∈A kalama b∈B kblbmb ˜ Qa kalamagn(Rab)Tklm(Rab)˜ Qb,xy kblbmb, (5) where k=ka+kb,l=la+lb,m=ma+mb,˜ Qa kalama are the cartesian multipoles of atom ain H2O expressed in the laboratory frame, ˜ Qb,xy kblbmbare the cartesian multipoles of atom bin OH for the matrix element Vxy expressed in the laboratory frame and Tklm(Rab) is the corresponding Coulomb interaction tensor.24 The symmetry constraint regarding the rotational invariance around the OH molecular axis is automatically fulfilled by this expression. This follows from the definition of Qb,xy kblbmb=x|xkbylbzmb|yin the molecular frame, which possesses necessarily the symmetry of the OH molecule. The function gn(R), with n=k+l+m, is a spherically symmetric damping introduced to avoid the unphysical divergence of expression (5) when Rapproaches 0. We use gn(R)=(1 −e−γR2)n.18 This function produces a shorter ranged damping than the popular Tang-Toennies damping. The H atoms have only a monopole here, while the O atoms were given a monopole, dipole, and quadrupole. The dipole and quadrupole of each molecule were calculated at the same level of the reference points. There is of course some freedom in the choice of the multipole parameters, provided that the total multipole of each molecule is reproduced. We constrained the oxygen charge of the OH molecule to be qB=−0.36. Since the dipole of the molecule is dOH =−0.65 a.u., the O atom carries a small dipole QO, 001 =−0.0124 a.u. along the symmetry axis of the molecule. The oxygen charge in the H2O molecule qA=−0.7622 was fitted on a set of long range data obtained for an oxygen-oxygen distance ROO >10.0 a.u. This allows to mimic approximately the small contribution of the missing high-order multipole not included in the fit. The values of the multipoles are given in Table II and in the supplementary material.18 The induction interaction is limited to dipole polarizability for the neutral molecules considered here. This interaction is not pairwise additive, and requires to compute the carte- TABLE II. Atomic cartesian multipoles used for the parameterization of H2O-OH interaction. For OH, Qxx ij k is given. The other cases are deduced by symmetry: Qyy ij k =Qxx jik and Qxy 110 =Qyx 110 =1 2Qxx 200 −Qxx 020. Only nonvanishing multipoles are given. The values Q∗ ij k are the actual molecular multipoles resulting from our parameterization. Molecule Center i j k Qijk Q∗ ij k H2OO000−0.7622 O001−0.1287 −0.72 O200−5.6300 −5.63 O020−4.7900 −4.00 O002−5.1700 −4.42 H 0 0 0 0.3811 OH O 0 0 0 −0.3600 O001+0.0124 −0.65 O200−4.0000 −4.00 O020−4.9800 −4.98 O002−4.2400 −3.17 H000+0.3600 sian components of electric field Eklm experienced by each molecule to form the induction matrix W. For OH there is one polarizability tensor αB,xy for each matrix element Wxy. The latter reads Wxy =− 1 2 klm  klm EklmαB,xy klm,klmEklm.(6) For H2O, there is one single tensor αA, but there is one contribution to the electric field from each (xy) arrangement and we have Wxy =− 1 2 klm  klm s=x,y Exs klmαA klm,klmEsy klm.(7) The electric field is obtained simply as the sum of the multipole contributions from all the surrounding molecules. With the above definitions, the electric field generated at Ra by the atom blocated at Rbreads Exy klm(Ra)= b klm gn(Rab)Tk+k,l+l,m+m(Rab)˜ Qxy klm, (8) where we have introduced the same damping function gn(Rab) as for the electrostatic contribution to avoid the short range unphysical divergence. The values of the cutoff were chosen so that the polarization is cut moderately at equilibrium distance, which ensures a better transferability to condensed phase. The value of the polarizability tensor were taken from Refs. 19 and 25. They are given in cartesian form in Table III. The dispersion interaction is limited to its leading isotropic contribution to R−6and R−8. The latter was introduced because its coefficient C8is relatively large with respect to C6.19 For the same reason, we also introduced dipolar C7 terms. There is no higher multipole moments included in our dispersion term, and therefore Vxy =0 and Vxx =Vyy.The complete expression for dispersion interaction reads Vxx =g6(R)C6 R6+g7(R)C7Aˆ R·uA R7 +g7(R)C7Bˆ R·uB R7+g8(R)C8 R8,(9) where Ris taken as the oxygen-oxygen distance, ˆ Ris the unit vector R/R, and uAand uBare the unit vectors oriented along the dipole of the molecules A and B, respectively. The values TABLE III. Cartesian polarizability tensors used for the parameterization of H2O-OH interaction. For OH, αyy ij k,ij k =αxx jik,jik and αyx ij k,i jk=αxy jik,jik. Only non-vanishing components are given. Mol. (st) i j k ijkαst ij k,i jk H2O1001009.30 010010 10.29 0 0 1 0 0 1 9.72 OH (xx) 1 0 0 1 0 0 6.37 0 1 0 0 1 0 7.75 0 0 1 0 0 1 8.75 OH (xy) 1 0 0 0 1 0 −1.05 010100 −0.13 Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44 164313-5 Galbis, Giglio, and Gervais J. Chem. Phys. 139, 164313 (2013) TABLE IV. Dispersion coefficient used in the fit. C6C7AC7BC8 −40 +50 +50 −1200 of the coefficient were computed from the dynamic polarizabilities given in the works of Meyer and Spelsberg.19,25 They are given in Table IV. The short range part of the interaction potential is dominated by hydrogen bonding. It is far from isotropic and for the sake of accuracy, we represent it as a sum of bond functions, which depends on the orientation of the bond with respect to the body-fixed frame of the molecules A and B. The matrix element expansion of Eqs. (2)–(4) are obtained from Vmdefined as follows: Vm= a∈A lama b∈B lb k calamablbmkfabk(R)Ylama(A)Ylbm(B) (10) with m=0, 2 and Ylm are real spherical harmonics. For the radial functions, we use a set of gaussian functions fabk(R)=e−kγabR2. The complete parameterization is given in the supplementary material.18 The choice of index for the spherical harmonics depends on the symmetry of the molecule. In the following, the symmetry axis of the molecule is aligned along the zcartesian axis. For hydroxyl (B), we can have only m=0, 2 and it was sufficient to limit lb≤2. For the water molecule (A), the C2vsymmetry constrains ma=0, 2 for the bond linked to the O atom. For the bond linked to the H atom, the parity of Ylamashould be such that the combination of the two bonds in the water molecule respects the C2v symmetry and thus ma=−1, 0, 2. It was sufficient to limit la ≤3 for the water molecule. The strong repulsion experienced by two very close nuclei is far above the energy of interest here and we simply represent this interaction by one repulsive function CABe−γABR/R for each pair of atoms. The latter interaction avoids in particular unphysical deep wells in the PES at very short distance. The values of the parameters are given in the supplementary material.18 It is important to note that, despite their apparent complexity, the computation of the matrix element for electrostatic, dispersion, and short range contribution reduces actually to a sum of pairwise interactions between the atoms of molecule A and those of molecule B. IV. FITTING PROCEDURE Before fitting the PES, a large set made of approximately 20 000 geometries was generated and the corresponding adiabatic energies calculated. In a first step, the long range interaction potential corresponding to distances ROO ≥10 atomic units was fitted by adjusting the oxygen charge in water and one exponentially screened function associated with 5 spherical harmonics on HAOBbonds, according to expression (10). In a second step, we retained the geometries for which the energy of the first PES is located in the range [−0.01; +0.01] a.u. around the asymptotic energy of the dissociated complex. FIG. 2. Energy difference between the fitted and reference PES as a function of the computed reference energy. For the second PES, some of the corresponding energies are much higher, in particular for short ROO. We retain thus approximately 16 000 geometries, associated each with one energy for each PES. The points were fitted by a standard leastsquare fit procedure with equal weights to ensure a correct balance between all the points retained for the fit. The global accuracy of the fit is illustrated in Figure 2, where we plotted the energy difference δ=Efit −Eref,asa function of the reference energies Eref. The standard deviation is better than 10−4a.u. (i.e., approximately 20 cm−1). Most of the reference points lay between the two lines δ±2 ×10−4a.u. The dispersion is slightly larger for the second PES for the highest energies. Such an accuracy is reached with a limited set of 242 parameters in the present case. This is rather modest regarding the complexity of the coupled PES and the relatively broad range of energy retained for the fit. Beside its intrinsic quality, the interest of a matrix formulation to represent the surface can be observed near the PES crossings or avoided crossings. In Figure 3, we have plotted a few cuts of the PES for a planar H2O-OH complex. Each PES is associated with a given symmetry Aor A of the Cssymmetry group and can thus cross each other. The Asymmetry corresponds to the singly occupied π-orbital of the hydroxyl in the symmetry plane, while the A symmetry corresponds to the singly occupied orbital perpendicular to this plane. The relative orientation of the molecules in the plane is given by θAand =θA−θB, where θAand θBare, respectively, the direction of ROO and the direction of the OH radical axis with respect to the symmetry axis of the water molecule (see Ref. 18). Thus, =0, corresponds to the hydroxyl aligned along the O-O direction and the hydrogen pointing outward. The scan is taken at ROO =5.5 a.u., which corresponds approximately to the equilibrium distance of the main minimum. As it can be observed in Figure 3, the fit is remarkably close to the control data, though the latter points were not included in the fitting procedure. In particular, the PES crossings are faithfully reproduced by our matrix potential approach, while they would be terribly difficult to obtain with a scalar fit of the adiabatic PES. Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44 164313-6 Galbis, Giglio, and Gervais J. Chem. Phys. 139, 164313 (2013) FIG. 3. Potential energy of the H2O-OH complex as a function of for some selected θAand ROO =5.5 a.u. From bottom to top: θA=0◦,θA= 38◦,θA=76◦,θA=104◦. The points are calculated control data and the lines correspond to the fit. V. PES ANALYSIS The impact of the PES crossings is not negligible as they are located around the secondary minimum, which controls the formation of hydrogen bond when OH acts as an acceptor and thus controls the motion of the atoms in the corresponding basin. A contour plot of both PES is shown in Figure 4for ROO =5.8 a.u. We can observe a series of minima in this PES cut. For both PES, the deepest minimum, on the right-bottom corner (=180◦, θA=0◦) comes from the saddle point between the two symmetric main minima associated with OH acting as a donor. For the ground state PES, the second minimum on the upper right part of the figure (≈80◦,θA≈110◦) comes from the secondary minimum associated with H2O acting as a donor. The figure shows clearly that this minimum is extremely shallow. For the second PES, the trace of the secondary minimum (≈20◦,θA≈110◦) is even more shallow in this plane. For this minimum, the equilibrium O-O distance is actually slightly larger than the value taken for this scan. The gap along the minimum energy path connecting the secondary minimum of a given PES to its associated main minimum is less than 5×10−4a.u. (100 cm−1) for the lowest PES and slightly larger for the second PES. For the ground state PES, we can observe another shallow minimum on the upper left side of Figure 4(≈−60◦,θA≈120◦). It is linked to the basin of the secondary minimum by rotation of the water molecule around the OH direction at constant θA. The main penalty comes from the multipole interaction, which is less favorable for this orientation. This seeming minimum is connected, without energy barrier, to the secondary minimum on the upper right side of the figure. The secondary minimum is thus characterized by two directions of easy motion along which the deformation of the structure would be controlled essentially by steric effects in condensed phase. Figure 4shows also the surface crossings, which can be followed by looking at the singular points along the isoenergy curves in the PES cut. At these points, an iso-energy line on the lowest PES continues by a line on the second PES and vice versa. One such crossings is visible around ≈20◦ Δ(degrees) -150 -100 -50 0 50 100 150 0 20 40 60 80 100 120 140 160 180 θA (degrees) FIG. 4. Potential energy contour of the H2O-OH complex as a function of and θAfor ROO =5.8 a.u. The red lines correspond to the lower PES and the blue lines to the upper PES. The iso-energy lines are separated by 2.5 ×10−3a.u. Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44 164313-7 Galbis, Giglio, and Gervais J. Chem. Phys. 139, 164313 (2013) for θA=20◦–100◦. Another one can be observed on the right hand side of the figure from ≈150◦and θA=20◦up to  ≈80◦and θA=70◦. A third crossing can also be observed on the left hand side of the figure for −40◦to −70◦and θA =30◦–60◦. At these crossings, the orientation of the hydroxyl quadrupole changes and thus influences the course of the dynamics in MD simulation. VI. TRANSFERABILITY To investigate the hydroxyl radical in solution, the H2O- OH intermolecular potential needs to be combined with a model for water molecule interaction. Our fit is combined here with the mobile charge density in harmonic oscillator (MCDHO) model.5Many alternative models taking explicitly the polarization into account and providing accurate radial distribution functions21,26 might be used instead of the MCDHO model. The polarization in the MCDHO model is based on an electric charge harmonically bound to the O atom of the water molecule and interacting by means of screened Coulomb interaction with the other charges of the MCDHO model. We modify slightly our model to make it compliant with the MCDHO model as follows. We simply replace the polarization of the water molecule (Eq. (7)) by the interaction of the hydroxyl multipoles with the mobile charge. For the sake of simplicity, we neglect the small difference in the water molecule polarization by the two OH diabatic states by taking the trace of the electric field tensor defined in Eq. (8).The polarization of the water molecule reduces thus to a diagonal matrix equivalent to the scalar potential energy used for a single PES model. The potential energy induced by polarization of a water molecule by the hydroxyl is simply the difference UB(qA)−UB(q0 A) where qAis the equilibrium position of the mobile charge associated with the water molecule A in solution, q0 Aits equilibrium position for the free molecule A at the same geometry, and UBthe electrostatic interaction energy of the mobile charge with the hydroxyl multipoles. The increase of internal energy by polarization comes simply from the interaction of the mobile charge with the atomic charge of the water molecule. The polarization of H2O becomes therefore the polarizability of the MCDHO model. In the simulation, the positions of all the mobile charges are obtained by minimization of the energy with respect to their positions. The corresponding self-consistent equations are solved iteratively for each move in the MC simulation. The modification of the energy functional and water polarizability with respect to our original fit functional leads to a small energy difference of the order of 2 ×10−4a.u. (roughly 40 cm−1) at the geometry of the main minimum. This difference originates partially from the difference between the water molecule polarizability tensors of the MCDHO model and the one used in our fit. To have the possibility to study the OH deformation in solution, the hydroxyl molecule is flexible in our simulation. Its internal energy is taken from a fit to the internal OH energy,27 which completes the representation of the two lowest PES of a system made of several water molecules. The modification of the OH length is actually very small because of the large difference in energy scale between the internal and external degrees of freedom. The resulting small perturbation of the multipoles is neglected here. We perform a transferability test by analyzing the (H2O)2OH trimer for various geometries. We take 100 geometries made of the two closest neighbors of OH from MC sampling of the liquid as described in Sec. VII.Wealsotake 50 configurations of the isolated trimer from MD runs of 10 ps each and starting at the equilibrium geometry with a set of velocities sampled from a Maxwell-Boltzmann distribution at a temperature T=300 K. Finally, we perform a scan by moving the OH along the direction uddefined by the positions of the hydroxyl oxygen and the middle of the two water oxygens at equilibrium. A common offset energy corresponding to the asymptotic position reached along this direction is subtracted to all the computed energies. We present in Figure 5, a plot of the correlation between the values of the fit and the corresponding ab initio energies computed with the MRCI method and ANO-S basis set described in Sec. II. The observed correlation is quite good for both states. We observe a small systematic deviation from perfect correlation. The fit underestimates the depth of the basin explored around the equilibrium geometry as can be observed from the points deduced from the MD runs and the scan along ud. The repulsion along this direction is largely overestimated for the shortest distances between the OH and the two water molecules. In this range, the PES is very steep and this large energy deviation corresponds actually to a rather small position shift, of the order of 0.02 a.u., which has thus a minor influence on the radial distribution function presented in Sec. VII. Moreover, this repulsive wall is not explored during MC simulation, as it can be deduced from the corresponding points plotted in Figure 5. The average dispersion around the perfect correlation line is of the order of 10−3a.u., which is comparable to the dispersion reported for the authors of the MCDHO model.5The FIG. 5. Correlation plot between the model and the corresponding ab initio energies for a selection of geometries for (H2O)2OH cluster. Filled symbols: state 1; hollow symbols: state 2. Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44 164313-8Galbis, Giglio, and Gervais J. Chem. Phys. 139, 164313 (2013) three-body interaction is not excellent for this model according to its authors.5A part of the observed deviation observed in Figure 5comes from the water-water interaction, and the use of alternative model, with a more consistent treatment of the polarization to represent this interaction, might produce a better result. Finally, it is interesting to note that the energy difference between the two adiabatic PES is quite well reproduced by the fit. This indicates that most of the error comes from the trace of the model and that the proposed diabatic approach has the proper functional form to represent the interaction of OH with water. VII. STRUCTURAL PROPERTIES OF OH RADICAL IN SOLUTION We investigate the properties of OH in liquid water by means of MC sampling with the Metropolis algorithm in the NVT ensemble. The density of the liquid is thus held constant at 1.0 g/cm3and the temperature is T=300 K. The liquid is represented by a cubic box containing 125 molecules with periodic boundary conditions and we used the minimum image convention. We sample the molecular moves by combining a translation of the molecule center of mass and a rotation around one of the axis chosen randomly. The largest amplitude for translation and rotation were 0.3 a.u. and 0.3 radians, respectively. The internal moves are sampled by combining a translation of each atom of the molecule. Since this motion is energetically much more expensive, the largest amplitude was limited to 0.02 a.u. As usual for MC simulation, the sampling amplitude was scaled in the course of the simulation to ensure a ratio of 50% between acceptance and trial. The results were obtained by sampling 2 ×108moves after a relaxation run of 107moves. Since the OH radical has two degenerate electronic states, it is apriorinecessary to consider both PES to sample the configuration of the liquid. Generally speaking, it is important to use a diabatic representation as we did, so that the contributions of the water molecules add coherently with each other. If only the adiabatic energy were added, the energy difference between the two PES would be unphysically too large. Considering the isolated H2O-OH complex, the energy difference between the two PES varies significantly with the geometry. In the area of the main minimum, it is of the order of 10−3a.u. (roughly 200 cm−1) and increases when the oxygen-oxygen distance diminishes. However, in a large portion of the configuration space, where OH acts as an acceptor, the energy difference between the two PES is significantly larger. For the configurations relevant for MC sampling in solution, it is likely to find at least one of the water molecules in such a position and the energy difference is of the order 8×10−3a.u. (approximately 1600 cm−1). This is roughly 8 times larger than kBTat T=300 K. In this case, the relative weight of the second PES is less than 10−4and it plays thus only a minor role. The oxygen-oxygen radial distribution function, gOO(R), around the OH radical is shown on top of Figure 6. It is quite similar to its counterpart for pure water. We observe a broadening and a small shift of the first peak position at 5.4 instead of 5.3 a.u. with the MCDHO model. The broadening FIG. 6. Radial distribution functions for OH in solution and for pure water. Continuous line: OH radical solvated in water. Dotted-dashed line: pure water. comes from the difference of the two minima in the H2O-OH PES with respect to the corresponding main minimum of the H2O–H2O PES. This peak does not present any sharp structure as reported in previous study based on a simplified PES.6 The position of the first minimum is 7.7 a.u. and we found approximately 9 oxygen atoms inside this first coordination shell. However, all these molecules are not hydrogen bonded to the OH radical. The hydrogen-oxygen distribution function, gHBOA(R), is depicted in the middle panel of Figure 6. It corresponds to the distribution of water (A) oxygen with respect to the hydroxyl (B) hydrogen. For R<8.0 a.u., the distribution is significantly different from the pure water distribution, while both distributions become identical for larger distance. We also notice a sharp contrast with the results reported by Cabral de Couto and co-workers,6who observed a well-defined sharp peak at short distance (R=3.1 a.u.). In our case the first peak is located around R=3.4 a.u., at a slightly shorter distance than for pure water. As noted previously for water,5the role of polarization is essential to ensure a correct transferability of the fitted PES from isolated bimolecular complex to condensed phase. A similar simulation neglecting the polarization produces a donor hydrogen bond centered around R=3.7 a.u. instead. The first peak is higher than its counterpart for pure water. This is mainly due to the presence of the second peak in the radial distribution function. The latter is centered around R=6.5 a.u. but it is significantly broader than its counterpart for pure water. The broadness reflects the shape of the H2O-OH PES around the secondary minimum, in which the Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Fri, 30 Sep 2016 13:47:44