Analytical free energy gradient for the molecular Ornstein-Zernike self-consistent-field method

An analytical free energy gradient for the molecular Ornstein-Zernike self-consistent-field (MOZ-SCF) method is presented. MOZ-SCF theory is one of the theories to considering the solvent effects on the solute electronic structure in solution. [Yoshida N. et al., J. Chem. Phys., 2000, 113, 4974] Molecular geometries of water, formaldehyde, acetonitrile and acetone in water are optimized by analytical energy gradient formula. The results are compared with those from the polarizable continuum model (PCM), the reference interaction site model (RISM)-SCF and the three dimensional (3D) RISM-SCF.


Introduction
The structure of molecule changes, sometimes drastically, upon solvation.Such a change plays an important role in determining the physical properties of the molecule, chemical reaction and biological processes [1].The optimal molecular geometry in solution has the lowest free energy in any other geometries.Therefore, in order to efficiently find the optimal geometry, free energy derivative with respect to the coordinate of atoms, called free energy gradient, is required.Although the free energy gradient can also be evaluated by numerical differentiation, such a procedure is not profitable in accuracy and computational time compared to the analytical method.
Considering these requirements, several theoretical methods have been proposed for analytical expressions of the free energy gradient.Such theoretical methods describe an electronic structure in a solution by combining the quantum chemical treatments for solute electronic structures with the classical statistical mechanical treatments for solvent distributions around the solute.One of the most popular theoretical models is known as a dielectric continuum model [2][3][4][5][6][7].Since the solvent is regarded as a continuum media in this model, the molecularity of solvent molecule could not be considered.For example, it is incapable of describing a local solute-solvent interaction such as hydrogen bonding.
In order to maintain the molecular aspects of solvents, a theoretical model, referred to as the reference interaction site model self-consistent-field (RISM-SCF) method, has been proposed by employing the RISM integral equation theory in the statistical mechanics of molecular liquids to obtain the solvent distribution around a solute [8][9][10][11].An analytical expression of the free energy gradients in the RISM-SCF formalism was proposed by Sato et al. [12,13].They have reformulated the RISM-SCF equations with the use of the free-energy expression introduced by Singer and Chandler [14].It has been shown that the derivation automatically provides variational conditions for multiconfigurational self-consistent-field (MCSCF) wave functions as well as the Hartree-Fock (HF) wave functions.The expression for the first derivatives of free energy with respect to solute nuclear coordinates has been derived by using the variational conditions.More advanced threedimensional (3D) RISM theory is also implemented to ab initio electronic structure theory by Sato, Kovalenko, and Hirata [15].That method, called 3D-RISM-SCF, has been successfully applied to N.Yoshida fully anisotropic solute molecule, even biomolecule, and the analytical expression of the free energy gradient in 3D-RISM-SCF formalism has been proposed [16].However, both these methods employ the site-site Ornstein-Zernike (OZ) equation for describing the solvent-solvent correlation, which provides only spherically averaged correlation functions depending only on the site-site distance, and is known to give trivial values for the solvent dielectric constant.
Most recently, the molecular Ornstein-Zernike self-consistent-field (MOZ-SCF) method has been proposed in our previous paper [17] employing the MOZ integral equation theory.The MOZ theory allows us to treat the orientation dependence of intermolecular interactions through the rotational invariant expansions [18][19][20][21].The MOZ theory turns out to be efficient in reproducing the thermodynamic, dielectric, and structural properties obtained from the simulations for aprotic solvents, though there is a controversy as to the accuracy of the HNC approximation for the dielectric constant of protic solvents [22,23].In the previous paper, the multidimensional integral method proposed by Lado et al., which is applicable of treating the anisotropic molecule, was employed to solve the HNC closure [24].The quantum mechanical effects, i.e., the exchange repulsion and charge transfer interactions, have also been incorporated in calculation of the solute-solvent interactions.Since our approach is based on a molecular model of solvent, the short-range exchange repulsion and charge transfer interaction between solute and solvent at a molecular level were explicitly considered by introducing an effective potential located on solvent molecule.
In this paper, the analytical expression is proposed based on the MOZ-SCF model for the free energy gradient with respect to solute nuclear coordinates.It is possible that the MOZ-SCF model has been derived based on the variations of the free energy functional, similarly to RISM-SCF.
The organization of this article is as follows.In the section 2, the MOZ-SCF formalism and the analytical energy gradient are described as well as the computational details.The section 3 shows the results of calculations.The optimized geometries of water, formaldehyde, acetonitrile and acetone in water are compared with the results of RISM-SCF, 3D-RISM and polarizable continuum model (PCM).The concluding remarks are given in the section 4.

Theoretical method 2.1. MOZ-SCF formalism
Since the details of the MOZ-SCF method have been presented in our previous paper [17], only the outline of the theory is described here.
The total Helmholz free energy A of the system is given as the sum of solute electronic energy and the excess chemical potential coming from the solute-solvent interaction, Here H 0 is Fock operator in gas phase, ψ is the wave function of solute, and ∆µ is the excess chemical potential which is the functional of the solute-solvent total correlation function h(ω 1 ω 2 R 12 ) and the direct correlation function c(ω 1 ω 2 R 12 ) as well as the solute-solvent interaction energy u(ω 1 ω 2 R 12 ) [14]; where (12) = (ω 1 ω 2 R 12 ).ω 1 and ω 2 are the Euler angles defining the orientations of the solute and a solvent molecule, and R 12 is the vector connecting the origin of the solute coordinate and the center of mass of a solvent molecule, respectively.Here . . .ω1ω2 means the integration over the orientation angles, ω 1 and ω 2 .The hypernetted chain (HNC) closure is assumed for relating the correlation functions and the interaction energy, where g is the radial distribution function, h + 1, and β = 1/k B T with k B and T being the Boltzmann constant and temperature.
In the MOZ-SCF approach, the correlation functions as well as the interaction potential are expanded in terms of the basis set of rotational invariants; where R12 denotes the orientation of R 12 .The rotational invariants have the usual definition, where R m µ µ (ω) is a Wigner generalized spherical harmonic, [25] the brackets (• • •) is a 3-j symbol, and f mnl is conveniently chosen to be [(2m + 1)(2n The free energy A is regarded as a functional of the correlation functions, c(12), h( 12) and η( 12), as well as the one particle orbital and CI expansion coefficients.Imposing the constrains to the orthnormality of the configuration state functions and one particle orbitals, we can define the following Lagrange function: where S im = φ i | φ m , C I denotes the CI coefficient and im is the Lagrange multiplier.Variations with respect to the functions yield 12)+η( 12) δη( 12) ω1ω2 Here Cmn µν,χ , Hmn µν,χ and Xmn µν,χ are χ transform of the rotational invariant expansion coefficients, c mnl µν , h mnl µν and x mnl µν , respectively [18][19][20].x mnl µν is the coefficient of the solvent-solvent total correlation function.φ, h, g kl , γ ij and Γ ijkl are one particle orbital, one and two electron hamiltonian and the vector coupling coefficients.The variation of L with respect to the correlation functions and electronic wave function gives the MOZ integral equation, the HNC relation and the solvated Fock equation.The wave function of solute molecule is determined by and Here r, ρ and ûelec are the coordinate of the electron, the solvent number density and the solutesolvent pair interaction potential operator, respectively.g(ω 2 R 2 ) = g(ω 1 ω 2 R 12 ) ω1,R1=0 means that the orientation of solute molecule is fixed and the solute molecule is located at the origin of the coordinate.Note that the angle bracket means angler averaging of solvent orientation.The solute-solvent interaction potential operator is given as the sum of the short range and long range parts, ûelec (r; The short range potential represents the exchange repulsion between the electrons in a solute and electrons in a solvent molecule, which is given by where The radial part of the effective one electron operator is represented by the sum of simple Gaussian functions, where Y nν (r) is a spherical harmonic.The long range part represents the electrostatic interaction is given as follows: Here Q n ν is solvent multipole moment, and Qm,elec µ is related to solute multipole moment via in which Qm,elec , and R i is position of ith nuclei.The total solute-solvent interaction potential , u in equation ( 3), is the sum of the electronic, nuclei-solvent electrostatic and dispersion interaction; Here u Disp denotes the dispersion energy, for which the London static polarizability approximation is employed; [26] u Disp (12 where Ē12 = E 1 E 2 /(E 1 + E 2 ) with E 1 and E 2 being the ionization energies of solute and solvent, respectively.α mµ is an element of the polarizability tensor.In this study, the ionization energies of solute and solvent and polarizability tensor in London formula are assumed to be constant.Lennard-Jones (LJ) type potential can be employed instead of the sum of the effective shortrange potential and the London dispersion mentioned above.Lennard-Jones potential is written as where σ and are usual Lennard-Jones parameters and R ij is the distance between solute site i and solvent site j.Thus, when the Lennard-Jones potential is employed, the total solute-solvent interaction potential becomes In this case, only the electrostatic potential depends on the solute electronic structure.

Analytical energy gradient
In order to perform the geometry optimization for solvated molecule, the analytical energy gradient is introduced by a similar procedure described in reference [12].The first derivative of the free energy with respect to the nuclear coordinate of the solute molecule R a is written as + −1 − h( 12) + e −βu( 12)+η( 12) ∂η( 12) ∂R a ω1ω2 where the notation of symbols is the same as in reference [17].u core (12) is the electrostatic interaction potential between solute nucleus multipole moment, Q m µ , and solvent multipole moment, Q n ν .Using the MOZ-SCF variational conditions, equation (7), it easy to show that the free-energy gradient is reduced to where the MO integrals H a ij , (ij|kl) a and S a ij are calculated by transforming the derivatives of the corresponding atomic orbital integrals.The derivatives of u core (12) and u elec,a ij (12) are easily calculated in the same way as the other derivatives.In this study, since the polarizability tensor and the ionization energies of solute molecule are assumed to be constants, the derivatives of the dispersion energy are zero.
If Lennard-Jones potential is employed instead of the sum of the effective short range potential and London formula, the derivative of Lennard-Jones potential should be added to equation (21) instead of the effective short-range contribution to the free-energy gradient,

Results and discussion
The analytical free-energy gradient of MOZ-SCF was applied to a geometry optimization of H 2 O molecule in water solvent.Before performing the MOZ-SCF calculations, the solvent pair correlation function was obtained by solving the MOZ equation for pure water with the HNC approximation.The temperature is 300K and the number density is 0.033327 molecules Å−3 .The solvent parameters were taken from SPC model potential [27].The maximum order of truncation of solvent rotational invariant expansion is four, n max = 4.The order of quadrature employed solving the closure is 7.For all MOZ calculations, the same solvent pair correlation function was employed.In all numerical calculations of MOZ, the number of grid points was chosen to be 512 with the width of 0.07 Å.
In order to perform the solute-solvent calculation, MOZ code was incremented to the GAMESS (US) program package [28].The HF wave function was employed with the 6-31G* basis set [29].The rotational invariant expansion was truncated at m max = 4.The multipole moment of solute H 2 O was taken up to the octapole moment.
The MOZ-SCF calculations were performed with two different type potentials; one is the sum of the effective short-range potential, London dispersion and the multipole electrostatic interaction, and the other is the sum of the Lennard-Jones potential and the multipole electrostatic interaction.In the later case, Lennard-Jones parameters of solute molecule were employed in the same way as solvent one.
The geometry optimization was also examined using PCM, RISM-SCF and 3D-RISM-SCF method to compare the results with the MOZ-SCF.The PCM calculations were carried out using the option of default water.Both of RISM-SCF and 3D-RISM-SCF calculations were employed using the same Lennard-Jones parameters and basis functions with MOZ-SCF calculation.In the RISM-SCF calculations, the number of grid points was chosen to be 2048 with the width of 0.05 Å.On the other hand, 128 grids points with 0.5 Å were employed for 3D-RISM-SCF calculation.In table 2, the results of geometry optimization of H 2 O molecule in water solvent are shown.The results of MOZ-SCF calculation with effective short-range potential (EP) for geometrical parameter, O-H bond length and H-O-H angle, are close to gas phase value rather than the results of MOZ-SCF with LJ potential.On the other hand, both calculations show the similar value of dipole moment.As seen in figure 1, u 000 00 s, the lowest order terms of rotational invariant expansion coefficient of the interaction potential for both cases are very close.There are small differences which are depth of well and gradient of potential wall.The lowest order terms of rotational invariant expansion coefficient of the distribution function, g 000 00 , of two types of MOZ-SCF results are described in figure 2. g 000 00 denotes the radial distribution which depends only on the distance between solute molecular center and solvent molecular center.Although the positions of 1st and 2nd peaks of g 000 00 are quite the same, height and width of 1st peak for EP are lower and broader than LJ case.This difference is caused by the gradient of u 000 00 at the contact distance.Since EP employs the Gaussian function as a basis to reproduce the short-range interaction, u 000 00 of EP shows gentler slopes than LJ which employs 1/r 12 .However, these differences are quite trivial.Actually, these two calculations yield similar coordination numbers, 9.56 and 9.52 for EP and LJ at 4.0 Å, respectively.As seen in table 2, the results of optimized geometry of MOZ-SCF, RISM-SCF and 3D-RISM-SCF are very close.Although PCM results show a similar value of H-O-H angle, the distance between oxygen and hydrogen is close to gas phase value rather than other solvated structures.These results indicate that the optimized structure is strongly affected by the short-range potential rather than by the type of integral equation.On the other hand, the dipole moment shows the obvious difference.Therefore, the electronic structure of solute molecule is affected by the type of integral equation.
The geometry optimization of H 2 CO is also performed with C 2 v symmetry in aqueous solution.The wave functions employed here are the complete active space (CAS) SCF, which is constructed by distributing six electrons in the five active orbitals, i.e., the carbonyl π, π * , σ, and σ * orbitals and the oxygen nonbonding orbital.The basis set is the 6-31G*.Lennard-Jones parameters of solute H 2 CO molecule were taken from OPLS-AA parameter set [30].In table 3, the results of geometry optimization of H 2 CO molecule in aqueous solution are summarized.The results show the tendency similar to H 2 O case.RISM-SCF gives strongly polarized electronic structure and PCM underestimates the polarization.MOZ-SCF and 3D-RISM-SCF exert an intermediate effect on the solute electronic state between the RISM-SCF and PCM.These are comparable to our previous results [17].
The optimized geometries and dipole moments of acetonitrile are also given in table 3. The C 3 v symmetry is employed for the electronic structure calculations and C ∞v for the MOZ calculation.Lennard-Jones parameters of solute acetonitrile were taken from OPLS-UA parameter set [31].The dipole moments as well as the electronic structure show a tendency similar to the H 2 CO case.Although the optimized geometries are very close, MOZ-SCF gives a slightly different value for the distance between two carbons.It often occurs that double/triple bond is weakened by solvation effect.In fact, C=O is elongated in all the cases as shown in table 3.However, C≡N bond is not.Only in the MOZ-SCF it gives the decrease of the bond length.
The calculations of geometry optimization of acetone in aqueous solution were also performed with C 2 v symmetry.The HF wave function was employed with the 6-31G* basis set [29].Lennard-Jones parameters of solute acetone molecule were taken from OPLS-UA parameter set for MOZ-SCF [31].As seen in table 3, the structure evaluated by MOZ-SCF with m max = 4 shows a different feature from those by the RISM-SCF, 3D-RISM-SCF and PCM.It might be caused by the fact that the accuracy of the rotational invariant expansion is not sufficient.Since anisotropy of acetone molecule is relatively high, the rotational invariant expansion is short to reproduce the solute-solvent correlation.The MOZ-SCF calculation with m max = 5 was also performed to study the effects of truncation of the rotational invariant expansion.The C-CH 3 bond length approaches to RISM results rather than MOZ with m max = 4, though the angle CCC becomes narrower.On the other hand, the dipole moments, as well as the electronic structure, of both MOZ results show the similar tendency to H 2 O and H 2 CO case.The convergence of the electronic energy seems to be fast, namely, only lower order contributions are sufficient for the description of the electronic energy.However, much higher terms are necessary for the accurate energy gradient.

Conclusion
In the present paper, the analytical expressions for the gradient of free energy with respect to the solute nuclear coordinates were derived from the variational conditions for the MOZ-SCF method.The calculations of the geometry optimization of H 2 O and H 2 CO molecule in water solvent were carried out with the HF and CASSCF wave functions.The optimized geometry with the present method was compared with those predicted using the RISM-SCF, 3D-RISM-SCF and The difference of the effect of short-range potential between LJ and EP was considered.The effect on the solute electronic structure turned out to be very close.Thus, these calculations gave similar dipole moments.On the other hand, the structure of solute molecule was strongly affected by the short-range potential.The results of H 2 O and H 2 CO show a similar tendency.The MOZ-SCF gave comparable results to those from other methods.
The MOZ-SCF showed the different trends in the optimal geometry of acetonitrile and acetone.It might be caused by the fact that the accuracy of the rotational invariant expansion is insufficient.The higher order terms of rotational invariant expansion were required to evaluate the optimal geometry in comparison with the electronic structure evaluation for anisotropic molecule such as acetone.Unfortunately, MOZ calculation with a higher limit of truncation of rotational invariant expansion demands much computational time compared with RISM and 3D-RISM.
The implementations of RISM-SCF/MCSCF and its analytical energy gradient have been successfully used in order to study the reaction dynamics of polyatomic systems, geometry optimization of photo-excited molecules and so on [32][33][34][35].Since the present model is derived based on the same logic as RISM-SCF, such studies with the MOZ-SCF method can be performed, which are in progress in our laboratory.

Figure 1 .Figure 2 .
Figure 1.Comparison of the center-center interaction u 000 00 of models for water.

Table 1 .
The parameters C n ν,k and ζ n ν,k were least-squares fitted to the ab initio results of the sum of exchange repulsion and charge transfer energies for solvent dimer.The parameters C n ν,k and ζ n ν,k for water are summarized at table 1. Potential parameters of effective one electron operator for H2O given in atomic units.

Table 2 .
Geometrical parameter as a result of optimization of water in water.and b are MOZ results with the effective short-range potential and with the Lennard-Jones potential, respectively. a

Table 3 .
Geometrical parameter as a result of optimization of formaldehyde, acetonitrile and acetone in water.Values in the parentheses are the results with lower order truncation, mmax = 4. a