Energy partitioning on intermolecular interactions: ab initio Monte Carlo study of water dimer

Ab initio Monte Carlo computations were carried out on H2O dimer system. By introducing the energy partitioning scheme that we have developed recently, ab initio calculated H2O–H2O interaction can be analyzed from the viewpoint of atom-atom interaction. The electronic polarization caused by the interaction and its temperature dependence are also discussed. To our best knowledge, this is the first report on the thermal distribution of electronic distortion energy assigned to a molecule.


Introduction
In aqueous solution, a strong intermolecular interaction exists between water molecules.It is the hydrogen-bonding that characterizes the structure of the liquid.In order to understand the nature of this ubiquitous liquid, a huge number of molecular simulations have been performed so far.Most of them are based on empirically derived intermolecular potential energy functions, which is essentially the sum of Coulombic interaction and repulsive interaction.
Of course we know very well that these functions can reproduce various properties of liquid water as a result.But how can we trust these empirical functions from the viewpoint of the first principle?In other words, can we construct the intermolecular potential energy functions from the first principle without any empirical adjustments?The electronic structure of a molecule is considerably changed when the molecule interacts with other species.Even if its electronic wave function is obtained, it is numerical description of the molecule, and is often far from our chemical sense or chemical intuition.Thus, a bridge between the intuition and modern computational results described by wave functions is needed.
One of the attempts to meet such demands is to partition the total energy of a molecule as follows, where E I is energy contribution from atom I in the molecule (one-center part), and E IJ represents energy contribution from the interaction between two atoms, I and J (two-center part).The energy partitioning has been well recognized in semi-empirical method.But it is not the case for ab initio molecular theory since the treatment of many-body contribution is not trivial.The way to partition ab initio total energy has unlimited number of possibilities, and there is no unique way of defining E I and E IJ .Ichikawa et al. proposed a partitioning scheme which describes bond-dissociation limit correctly [1].Mayer reported somewhat different scheme [2].Nakai has proposed another type of partitioning scheme called bond-EDA [3,4].
Recently we have proposed another scheme [5].As mentioned above, there are innumerable possibilities to introduce energy partitioning scheme for the electronic energy, but it should be emphasized that our method provides physically reasonable behavior for a wide range of interaction potential energy, and the intermolecular interaction is treated equivalently with intramolecular interaction.We also found that the components of molecular interaction in aqueous solution closely resemble those obtained by widely used empirical intermolecular potential functions such as TIP3P.
In the present study, we combined this partitioning scheme with ab initio Monte Carlo technique.The water-water interaction is undoubtably the fundamental interaction in condensed phase, and the minimum-size molecular pairs having hydrogen-bonding.Since ab initio molecular orbital calculations should be repeated at each step of the Monte Carlo simulations, together with the partitioning analysis, the computational costs is no small matter even for this small system.The present computation can offer a new view of the intermolecular interaction based on the site-site description, completely from the first principle.It should be emphasized that the electronic distortion energy assigned to a molecule can be computed only by introducing energy partitioning scheme.To our best knowledge, this is the first report on the thermal distribution of such distortion energy.

Energy partitioning scheme
We shall start with the total energy of the system, where h ab is the matrix element of the one-electron hamiltonian over the basis orbitals {χ a }, P ab is the "density matrix" where f i ≡ 2 for closed shell system.Two electron integrals are defined as follows, Note that all the centers of the basis functions are located on any one of the atoms composing the molecule.Using the symmetry of the two electron integrals of the basis sets in real numbers, we finally reach the following expression, (ab|cd) P ab P cd − 1 2 (ab|cd) P ab P cd − 1 2 As shown in the previous study [5], the present expression provides a physically reasonable behavior for the energy components.It is also noted that the following reduction to one-center contribution of this expressions, comes down to the energy defined in EDA [6].
In the present study, we are interested in the intermolecular interaction.Based on the present partitioning procedure, we can define the interaction between a water (W1 consisting of three atoms, O1, H1 and H1 ), and the other water (W2 including O2, H2 and H2 ), as follows: The remains of the total energy express the locally defined energy of each molecule.
When the molecules interact with each other, the change in the energy of the molecule is evaluated with respect to that of an isolated water molecule, E W0 .
Note, however, that the label (W1 or W2) is arbitrary and there is no physical meaning to distinguish the water molecules in the actual Monte Carlo step.The following quantity is much more practically useful, By considering the total energy of the dimer system (E W1+W2 ), the interaction energy computed by the standard Hartree-Fock procedure (∆E) is related to these quantities.
It should be remembered that E X pol represents the deviation of electronic distribution in individual water molecules from their isolated situation, and does not include contribution from the interaction between them after this electronic change.E int expresses the interaction energy between the electronically altered molecules.Note that both of E pol and E int are not constant values but depend on the relative configuration of water.The distribution functions of these quantities are computed as described below.

Ab initio Monte Carlo computations
Standard metropolis algorithm was employed to generate the water-dimer configurations.At each step, the energy was computed by ab initio molecular orbital theory at the level of Hartree-Fock method.6-31G(d,p) basis sets were used throughout the study.The geometry of each water molecule is fixed at the optimized one (R OH =0.943 Å and ∠HOH = 105.968• ) since we are focusing on the intermolecular interaction.If the molecular geometry is flexible, the partitioned energy assigned to a specific atom can be attributed both to the change of the molecular geometry and to the change in the configuration between the molecules.All the computational procedures were carried out with program code GAMESS [7] modified by us.
Sampling was collected until the sufficient ensemble was obtained (∼ 2.0•10 5 ) at the temperature of 50 K, 100 K and 150 K, respectively.The intermolecular interaction becomes too weak at 200 K or higher.The entropic contribution overcomes the interaction between the molecules and they begin to dissociate each other.This phenomena can be easily controlled by introducing a biased sampling technique.However, in the present study, we are interested in behaviors at relatively low temperature, because of the connection to the standard ab initio molecular orbital calculations corresponding to 0 K.
The maximum of the obtained distribution is located around It is well known that the potential energy minimum of water dimer is this linear configuration.Meanwhile, the broadening of the distributions depends on the temperature.The distribution is strictly localized 2.7 Å< R OO < 3.2 Å and θ < 0.5 at 50 K, while it expands to 2.5 Å< R OO < 3.7 Å and θ < 1.0 at 150 K.As described in equation ( 8), there are four interacting O• • • H pairs in the dimer system.Å are assumed (see the first peak position in figure 2).In the similar manner, when the distance is 3.4 Å, corresponding to the second peak in figure 2, the interaction energy is estimated about −31 kcal/mol.At the minimum energy structure of water dimers, the distances of remaining three O• • • H pairs are about 3.4 Å (∼ -30 kcal/mol) and 3.9 Å (∼ −27 kcal/mol).While they are slightly lower in energy than the weaker interaction shown in figure 1, they are reasonably close [9].In summary, the interaction energy given by equation ( 6) is very similar to the classical estimation although the present partitioning scheme deals with the interaction described in the electronic wave function of the system.

The hydrogen-bonding and its atomic level partitioning
As the temperature raises, the distribution of the hydrogen-bonding interaction shown in figure 1 becomes broader and the peak position shifts to the right-hand side, say weaker interaction, which is consistent with the depression of the first peak in the RDF.On the contrary, the distributions of the weaker interaction energy (other three O• • • H pairs) do not show significant temperature dependence.At the low temperature of 50 K, the two peaks found both in EDF and in RDF are separated from each other, indicating that the hydrogen atom making hydrogen-bonding keeps the bonding through to the end.In other words, the configuration of the molecules is extremely localized near the potential energy minimum and the 'roles' of the atoms are rarely exchanged.As increasing the temperature, the two peaks in EDF and RDF begin to connect with each other, and the hydrogen atom making the hydrogen-bonding can be replaced with other atoms much more easily.
One should notice that this interaction energy of O• • • H pair (∼ −50 kcal/mol and ∼ −30 kcal/mol) is much stronger than standard value of hydrogen-bonding energy (∼ −5 kcal/mol).It is obvious that the another contribution is necessary to complete the hydrogen-boning.Figure 4 shows the interaction energies between the water molecules, ∆E and E int .∆E is the interaction energy between the molecules, which is commonly recognized within the framework of ab initio molecular orbital theory.As mentioned above, the energy is found around −5kcal/mol.Since the left-hand-side edge of this distribution corresponds to the bottom of the potential energy of the interaction, the minimum limit is absolutely fixed and the distribution becomes narrower as the temperature decreases.The width of the distributions are directly related to the width of the effective interaction potential well, in which all the relative orientations are integrated over.

Interaction energy between water molecules
E int is negatively greater by about 3kcal/mol than ∆E because the distortion energy of the electron clouds (E pol ) is not included in E int .
The width of the distribution is slightly wider than that of E int .Note that the fineness of the horizontal scale is much smaller than the above figures.The overall character of the distribution and its temperature dependence are similar both in ∆E and E int .E pol is plotted in figure 5, in which two separated peaks are clearly recognized at any temperature.One is located around −1 kcal/mol and the other is around +4 kcal/mol.It should be noted that this component is related to the electronic distortion (relaxation) with respect to that of isolated molecule, as well as to the number of electrons belonging to the molecule.If the number of the electrons is conserved with the same geometry, the distortion always raises the energy positively, due to the variational principle for the wave function.If the number of the electron is not conserved, it can contribute to either increase or decrease of the energy.The two peaks crossing the zero energy suggest that the electron is more or less interchanged between the molecules through the hydrogenbonding.Actually, the effective charge assigned to the proton donor water is about −0.03|e| (Mulliken or Löwdin population), and E W of this water (equation 9) is lower than that of isolated molecules.In other words, the left-hand side peak in the figure corresponds to the proton donor water while that in the right-hand side is the proton acceptor.At 50 K these peaks are completely separated from each other, indicating that the roles of the water molecules in the hydrogen-bonding is always fixed.However, as temperature increases, the two peaks begin to overlap, and the water molecules can act as both the donor or acceptor in the hydrogen-bonding.It is also noteworthy to point out that the midpoint of these centers is not located at 0.0 kcal/mol but at nearly 1.5 kcal/mol.Because of the variational principle for the electronic structure, the distortion energy of the overall dimer system should be positive.

Conclusions
In the present study, hydrogen-bonding, which is the representative intermolecular interaction in aqueous solution phase, is investigated from the viewpoint of locally defined energy.By using the energy partitioning scheme combined with the Monte Carlo method, intermolecular interaction, the electronic distortion of a water molecule and their distribution were discussed.We would like to emphasize that the combination used in this study enables us to compute the distribution of E int and E pol for the first time, while the standard ab initio calculated Monte Carlo simulation can only compute the distribution of ∆E.
We found that the hydrogen-bonding between water molecules seems to be understood essentially as the Coulombic interaction.In other words, this traditional view of the interaction can be obtained from the first principle bridged by the energy partitioning scheme.The present method may be a useful tool for analyzing and systematic-understanding of the empirical potential energy parameters.
3.0 Å of O• • • O distance (R OO ) [8] and around 0.0 radian of the smallest intermolecular O• • • O-H angle (θ) at any temperatures (not shown here).Namely, the oxygen atom and one of the four O-H bonds are in co-linear position.

Figure 1 .
Figure 1.Energy distribution function of the intermolecular O• • • H pairs.

Figure 2 .
Figure 2. Radial distribution function of the intermolecular O• • • H pairs.

Figure 1
shows the energy distribution function (EDF) of intermolecular O• • • H pairs.At 50 K, two peaks are found around −53 kcal/mol and −23 kcal/mol.In view of the radial distribution function (RDF) of the O• • • H pairs in this system (figure 2), the stronger interaction (−53 kcal/mol) can be attributed to the direct interacting hydrogen-bonding pair.It is noted that the classical Coulombic interaction is estimated to −51 kcal/mol when two point charges (−0.8 and +0.4; typical effective charge of oxygen and hydrogen atom in water, respectively) separated from 2.1

Figure 3 .
Figure 3. Energy distribution function of the intermolecular O• • • O pairs.

Figure 3
Figure 3 shows the EDF of intermolecular O• • • O pair.Strong repulsive interaction appears around 60 kcal/mol.This is again consistent with the estimation of the classical Coulombic energy (∼ +70 kcal/mol) when assuming the effective charges separated from 3.0 Å (O• • • O distance).Temperature dependence of the distribution is similar to the O• • • H pairs: the peak shifts to the left-hand side and its height is depressed as the temperature increases.Thus, it might be possible to draw the following understandings: the hydrogenbonding between water molecules can be essentially described by the classical Coulombic interaction.Relatively great repulsive interaction between O• • • O and great attractive interaction between O• • • H compensate each other, and the resultant energy of the hydrogen-boding is onedigit smaller, say about −5kcal/mol.This would be the reason why an empirically derived simple point-charge-type model works surprisingly well to describe various properties of water.

Figure 4 .
Figure 4. Distribution of Eint and ∆E of the water dimer.

Figure 5 .
Figure 5. Distribution of E pol of the water dimer.