Polyatomic molecules in condensed phase : bond order index and solvation energy studied by RISM-SCF theory

The reference interaction site model – self-consistent field (RISM-SCF) theory was applied to a series of polyatomic molecules in aqueous solution, including CH4, NH3, H2O and HF. The bond order index, which characterizes the chemical bond, and its change by the solvation were studied by comparison with a dielectric continuum model. The change in the electronic structure was also associated with the solvation free energy. We find that the bond order indexes are not sensitive to the solvation effect even if the charge assigned to the atom is considerably changed. The distortion in the electronic energy is in proportion to the change in the solvation energy, showing that the linear response regime is a good expression of the solvation process examined here, in spite of solving the non-linear RISM-SCF equation.


Introduction
It is well-known that the electronic structure of a molecule is significantly changed when the molecule is dissolved into solvent.In the last few decades, theory for the electronic structure of solvated molecules has attracted a lot of researchers' attention, and various types of combination methods have been developed.These can be roughly classified into three categories, namely, the dielectric continuum model such as polarizable continuum model (PCM), quantum mechanicalmolecular mechanics (QM-MM), and the method based on the integral equation theory for liquids.RISM-SCF [1,2], which is a representative of the third category, compiles two fundamental methods in theoretical chemistry: one is the reference interaction site model (RISM) [3,4], and the other is ab initio molecular orbital theory.The method determines the electronic structure and the statistical solvent distribution around the solute in a self-consistent manner.It is a remarkable advantage that the RISM-SCF method can provide information on the microscopic solvation structure based on the statistical mechanics.We would like to emphasize that the RISM-SCF theory can be regarded as an extension of RISM theory as well as that of ab initio molecular orbital theory.The RISM-SCF has been successfully applied to numerous molecular phenomena including chemical reactions, chemical equilibria, charge transfer processes and so on [5].
The aim of this review is to characterize the chemical bond in condensed phase, which is definitely one of the central subjects in physical chemistry.The solvation effect by all means modifies the character of a chemical bond in the solute molecule, and, at the same time, this modification in the electronic structure affects the structure of the solvent molecules surrounding the solute.For the system we choose a series of hydrogenated compounds of the second period elements in aqueous solution, which are the most fundamental and ubiquitous molecules in the universe.Since all of them have only 'single' bond and their electronic structures in the gas phase are well known, they are thought to be the most suitable for a systematic comparison for the preset purpose.

Computational method
The RISM-SCF method compiles ab initio electronic structure theory and statistical mechanical theory of molecular liquid.Total energy of the solvation system is defined as the sum of quantum chemical energy of the solute (E solute ) and solvation free energy (∆µ) [2].
Since the electronic structure of the solute and the solvation structure around the solute are determined in a self-consistent manner, the wave function of the solute molecule is bending from that in isolated state.The energy difference between the solute in isolated state (E isolated ) and that in solution phase (E solute ) is a quantity to measure the contribution of "solvation effects" to the electronic structure.
where |Ψ and |Ψ 0 are wave functions in solution and in gas phase, respectively.H is the standard hamiltonian for the electronic structure of the solute molecule.
where the variables have their standard meanings.Solvation free energy within the framework of the present theory (excess chemical potential derived with hyper-netted chain approximation) is given by, h and c are the total and direct correlation functions, respectively.β = 1/k B T , where k B and T are the Boltzmann constant and temperature, respectively.The Fock operator of the RISM-SCF theory (F solv ) contains a solute-solvent interaction, V , The interaction energy between the solute and solvent reaction field is then given by, where b α is a proper population operator for atom α in the solute molecule and V α is the electrostatic potential acting on this atom.
ρ is number of density of solvent, g αs is pair correlation function between solute α and solvent site s, and q s is effective charge assigned to s.In the RISM-SCF theory, the electronic structure of the solute determined by equation ( 5) and RISM equation are coupled through equation (7).By solving the equation in a self-consistent manner, many-body effect on the solvation process is naturally taken into account.
All atom-type interactions were employed, whose Lennard-Jones parameters were taken from literature and summarized in table 1.The SPC-like water [6] was used to describe the solvent.All the van der Waals interactions between solute and solvent were determined by means of the standard combination rule.The density of solvent water was assumed to be 0.033337 molecule/ Å3 at 298.15 K.The PCM computations were carried out with standard parameters installed in the GAMESS program package [7].The electronic wave functions were computed based on RHF method with 6-31G(d,p) basis set [8].Molecular structures (geometry) of the solute molecules (CH 4 , NH 3 , H 2 O and HF) used in the present study were optimized using the standard RHF method without solvation effect.Computations in aqueous solution phase were then carried out by RISM-SCF method and by PCM method at this geometry.This is because we wish to genuinely know the changes in electronic energy due to solvation effect.Note that the optimized molecular structures in solution phase obtained by RISM-SCF and PCM methods are different.The electronic relaxation ascribed to the solvent effects and that ascribed to the difference in molecular structures (between RISM-SCF and PCM) are mixed up if the structural relaxation of the solute molecule is permitted.By this treatment, contributions from structural difference in the electronic structure and the related energy can be excluded.Upon transferring a molecule from gas phase to solution phase, the electronic structure of the molecule is altered by the solvation effects.In other words, the character of the chemical bond in solvated molecule is somewhat different from that in the gas phase.Although all of the chemical bonds examined in this study are formally referred to as "single" bond, the bond, its strength and order can be characterized in different ways.There is no unique definition to the evaluation of the bond order, and several definitions have been proposed.The definition proposed by Mayer, B Mayer , is considered as the most popular one [9].

Chemical bonds in aqueous solution
where P is 'P-matrix' (density matrix) and S is overlap matrix.µ and ν are the index of atomic basis sets running over the atom A and B. By analysing the molecular energy, we have recently found out that the following quantity could be another definition of the bond order index [10].
µ∈A ν∈B It should be emphasized that these both definitions are related to the covalent character of the bonding.As the ionic character of the bond increases, the bond order index becomes smaller.
Figure 1 shows the bond order index evaluated in the series of polyatomic molecules.The horizontal axis indicates the index of the molecule in the gas phase.As expected, the index becomes smaller with the increase of the electronegativity: the bond order indexes in methane are close to one, while those in HF are found to be around 0.85.The two definitions give very similar results for all the molecules although the difference is slightly enhanced in polar bonds such as H-F and H-O.The solvation effect is expected to change the character of chemical bond.The vertical axis of the figure corresponds to the bond order index in aqueous solution.The indexes of C-H bond in the gas and in aqueous solution are of the similar value, suggesting that the solvation effect is almost negligible in this molecule.This tendency on the C-H bond does not depend on the definitions of the indexes.The situations are slightly different in other molecules and the indexes decrease when the molecule is dissolved.The differences from those in gas phase become clear as the polarity of the bond increases.In other words, the covalent character of the bonds tends to be weakened by solvation effect.The indexes decrease as the polarity of the bonds increases, but the magnitudes of the change are only a few hundredth of a percent.In this study, two solvation theories are examined, i.e., RISM-SCF and PCM.Both of them provide very similar results although the solvation effect evaluated by the latter method is relatively weaker.Table 2 lists the Mulliken charges assigned to the hydrogen atoms.It is interesting to note that the charge assigned to the atom considerably changes due to the solvation effect although the bond order indexes are not sensitive to the solvation.Both of the electronic structures computed by RISM-SCF and PCM are polarized compared to that from the standard gas-phase computation, but they are slightly different from each other in detail.The PCM results are intermediate between RISM-SCF and gas-phase results.This tendency can be attributed to the same discussion presented in our previous study of a chemical reaction system [11].A localized interaction between specific atoms, such as hydrogen bonding, is explicitly treated in RISM-SCF method, while it is not in PCM.This difference sometimes plays a crucial role in determining the electronic structure of the solvated molecule.We should like to emphasize that the bond order indexes always decrease in an aqueous solution, and this conclusion does not depend on the choice of the definition of the index.It can be also said that the two solvation theories, RISM-SCF and PCM, offer a similar physicochemical picture of the solvation although they are quantitatively different.

Solvation free energy
In the present study, the electronic structure of the solute molecule is different from that in the gas phase, and the electronic structure and the solute-solvent interaction are determined in a self-consisent manner.On the RISM-SCF procedure the electronic structure is set to the gas phase wave function (|Ψ 0 ) at the beginning, and it is altered by introducing mutual interaction between solute electronic structure and solvation structure (total correlation function).The converged wave function in solution phase (|Ψ ) is usually polarized compared to that in the gas phase.Both of the wave functions, |Ψ 0 and |Ψ , are respectively related to the effective charge assigned to atoms, q 0 and q.The electrostatic interaction energy, which is the dominant contribution in the interaction between solute and solvent, is computed as follows; V 0 and V are the reaction field acting on the solute molecule introduced in equation ( 7).Each of them is respectively related to the solvation free energy, ∆µ 0 and ∆µ, corresponding to the energy before/after the RISM-SCF iteration cycle.Here we should like to recall the readers' attention to the fact that V 0 and V are directly related to the solvent distribution (see equation ( 7)), and their difference is considered to be the relaxation of the surrounding media due to the solute-solvent coupling.At the same time, q 0 and q are related to the electronic structure of the solute molecule.
In the RISM-SCF cycle, these quantities are coupled to each other and the many body effect is naturally taken into account.In what follows, the details of the changes due to the coupling are discussed.Figure 3 illustrates the change of the solutesolvent interaction on the RISM-SCF cycle.At the beginning, the wave function of the solute molecule in the isolated state (|Ψ 0 ) is related to the effective charge (q 0 ).The solvent distribution, the reaction field (V 0 ) and the energy (A 0 and ∆µ 0 ) are consistent with them (shown in open circles).By applying the self-consistent procedure to the coupling between solute and solvent, all the quantities are changed and finally converged to q, |Ψ , V, A , ∆µ and so on (filled circles).As seen in the previous section, the change on methane is negligibly small.The electrostatic interaction energy and the solvation free energy virtually remain at the state corresponding to q 0 even after the RISM-SCF iteration is converged.The situation is quite different in other three molecules.For instance, water molecule gains -18 kcal/mol of electrostatic interaction energy and -2 kcal/mol of solvation free energy when the gas-phase electronic structure (|Ψ 0 or q 0 ) is assumed.After the RISM-SCF iteration, the electronic-structure change in H 2 O leads to the -28 kcal/mol of electrostatic interaction energy and -7 kcal/mol of solvation free energy due to the induced solute-solvent interaction.The order of the change among the three species is directly related to the order of polarizability, α(NH 3 ) > α(H 2 O) > α(HF).
Solvation process modelled in the present study is indeed very complicated: solvation structure (or total correlation function) is obtained after solving non-linear RISM equation.The electronic structure (or wave function) is determined by the self-consistent calculation of Hartree-Fock type equation.One might think that a simple explanation of the solvation process is impossible.However, the framework of the process is understandable from the viewpoint of linear response regime.
Within the present RISM-SCF framework, the total energy of the system is given by equations (1) and (2).A can be thus rewritten as follows, Here we follow the protocol by Gao who analysed QM-MM energy [12].The total interaction energy between solute and solvent is given by, The quantity can be further decomposed into the following three contributions.
where the first term (∆µ 0 ) is the solvation free energy between the 'unpolarized' solute and the solvent, the second term (E reorg ; see equation ( 2)) is the same as E dist defined in the Gao's paper and the third (∆∆µ) is a net gain in the solvation free energy between the polarizable solute and the bulk solvent, which is equivalent to E stab .It is noted that the original analysis by Gao uses the interaction energy instead of the solvation free energy used here.According to the classical linear response theory, the gain in the interaction energy appearing in the third term is twice the energy penalty for distorting the solute electronic structure (E reorg ).
In the dielectric continuum theory, the solvation free energy is one-half of the electrostatic interaction energy [13], and the deviation from the state corresponding to 'unpolarized' solute molecule is then given by, Combining with equation ( 15), the following relationship may be obtained.Figure 4 plots the change in solvation free energy (∆∆µ) and that in the electrostatic interaction energy (∆E el ) as functions of the electronic distortion energy (E reorg ).All the quantities here were computed using the RISM-SCF method.The solid and dotted lines respectively represent the relationships, where ∆E U−V is the variable along the vertical axis.Clearly, the linear response regime is the good expression to describe the relationship between the electronic distortion and the gains in the solvation.We should like to emphasize again that RISM, Hartree-Fock and RISM-SCF are very complicated non-linear equations.By solving these equations faithfully, many-body effects are taken into account.However, the results shown here suggest that the non-linearity is virtually negligible in these solvation processes.Even in the system containing transition-metal complex, where strong electrostatic interaction is considered to exist due to the great charge on the metal, the solvation process can be understood within the same regime [14].

Conclusions
In the present study, typical polyatomic molecules in aqueous solution were systematically studied using RISM-SCF method, which is the combination of two fundamental theories, i.e., RISM and ab intio molecular orbital method.We found out that the bond order indexes were not sensitive to the solvation effect even if the charge assigned to the atom was considerably changed.The distortion in the electronic energy was in proportion to the change in the solvation energy, showing that the linear response regime was a good expression of the solvation process examined here in spite of solving the non-linear RISM-SCF equation.
The RISM-SCF theory offers a wide range of information on solvation systems, including electronic structure of the solute molecule (dipole moment, bond order index, effective charge and so on) as well as quantities in statistical mechanics.

Figure 1 .
Figure 1.The bond order indexes in gas phase and in aqueous solution phase.

Figure 2 .
Figure 2. The change of the bond order indexes upon transferring from gas phase to aqueous solution.

Figure 3 .
Figure 3.The change of solute-solvent interaction before and after RISM-SCF procedure.

Figure 4 .
Figure 4.The electronic distortion and the change in solvation energies.See the text for details.

Table 2 .
Mulliken charge assigned on the hydrogen atoms (|e|)