19

The RISM-SCF/MCSCF Approach for Chemical Processes in Solutions

Fumio Hirata and Hirofumi Sato

Institute for Molecular Science, Okazaki National Research Institutes, Okazaki, Japan

Seiichiro Ten-no

Nagoya University, Nagoya, Japan

Shigeki Kato

Kyoto University, Kyoto, Japan

I.SOLVENT EFFECT ON CHEMICAL PROCESSES

Most chemical reactions in nature as well as in the laboratory take place in liquid solutions. Chemical reactions of molecules in living systems take place exclusively in the solution phase. One of the most important tactics in organic chemistry is to choose an appropriate solvent for their reactions or to get higher yields. The reaction yield and rate are controlled by the solvent through changes in the free energy difference between reactants and products and thus in the activation free energy. In living cells, chemical reactions are controlled by the solvent with an additional complication—conformational fluctuations in biomolecules. In all those reactions, the solvent effect manifests itself not only through solute– solvent interactions and solvent reorganizations but also through intramolecular processes that are always associated with changes in electronic structure. In this regard, theoretical investigations of chemical processes in solutions are inevitably coupled with studies of quantum and statistical mechanics.

One of the most straightforward realizations of such couplings is the ab initio molecular dynamics (MD) approach originated by Car and Parrinello (CP) and published by Laasoninen et al. [1]. The method, which consists of solving the Kohn–Sham density functional (DF) equation by simulated annealing, has been proven to be a capable theory in many applications. However, the treatment has serious limitations or difficulties in terms of system size and the handling of excited states. In a sense, the approaches based on the path integral technique share similar difficulties with the CP theory for chemical processes in solution. Combined with molecular simulations and the statistical mechanics of liquids, the method could have revealed many interesting physical aspects of quantum processes in solution. But, again, it is largely limited to a simple and/or small system such as a harmonic oscillator or a solvated electron. In this respect, it is highly desirable to exploit the molecular orbital (MO) theory for the electronic structure, which has been proven to be the most powerful tool for exploring molecular processes. This method does not have the serious limitations characteristic of the CP and path integral approaches.

417

418	Hirata et al.

However, theories that are based on a basis set expansion do have a serious limitation with respect to the number of electrons. Even if one considers the rapid development of computer technology, it will be virtually impossible to treat by the MO method a small system of a size typical of classical molecular simulation, say 1000 water molecules. A logical solution to such a problem would be to employ a hybrid approach in which a chemical species of interest is handled by quantum chemistry while the solvent is treated classically.

Roughly speaking, three types of hybrid approaches have been proposed depending on the level of coarse graining with respect to solvent coordinates. Methods based on continuum models smear all solvent coordinates and represent solvent characteristics with a single parameter, namely, a dielectric constant. Molecular simulations, on the other hand, take all solvent coordinates into account explicitly and sample many solvent configurations in order to realize meaningful statistics for the physical quantities of concern. The approach based on the statistical mechanics of molecular liquids coarse-grains solvent coordinates in a level of the (density) pair correlation functions. In what follows, we briefly outline these three methods.

A. Continuum Model

The continuum model, in which solvent is regarded as a continuum dielectric, has been used to study solvent effects for a long time [2,3]. Because the electrostatic interaction in a polar system dominates over other forces such as van der Waals interactions, solvation energies can be approximated by a reaction field due to polarization of the dielectric continuum as solvent. Other contributions such as dispersion interactions, which must be explicitly considered for nonpolar solvent systems, have usually been treated with empirical quantity such as macroscopic surface tension of solvent.

A variety of methodologies have been implemented for the reaction field. The basic equation for the dielectric continuum model is the Poisson–Laplace equation, by which the electrostatic field in a cavity with an arbitrary shape and size is calculated, although some methods do not satisfy the equation. Because the solute’s electronic structure and the reaction field depend on each other, a nonlinear equation (modified Schro¨dinger equation) has to be solved in an iterative manner. In practice this is achieved by modifying the electronic Hamiltonian or Fock operator, which is defined through the shape and size of the cavity and the description of the solute’s electronic distribution. If one takes a dipole moment approximation for the solute’s electronic distribution and a spherical cavity (Onsager’s reaction field), the interaction can be derived rather easily and an analytical expression of the Fock operator is obtained. However, such an expression is not feasible for an arbitrary electronic distribution in an arbitrary cavity fitted to the molecular shape. In this case the Fock operator is very complicated and has to be prepared by a numerical procedure.

Numerous attempts have been made to develop hybrid methodologies along these lines. An obvious advantage of the method is its handiness, while its disadvantage is an artifact introduced at the boundary between the solute and solvent. You may obtain agreement between experiments and theory as close as you desire by introducing many adjustable parameters associated with the boundary conditions. However, the more adjustable parameters are introduced, the more the physical significance of the parameter is obscured.

B. Simulations

Molecular simulation techniques, namely Monte Carlo and molecular dynamics methods, in which the liquid is regarded as an assembly of interacting particles, are the most popular

RISM-SCF/MCSCF for Processes in Solutions

419

approach. Various properties related to the macroscopic and microscopic quantities are computed by these methods.

A typical hybrid approach is the QM/MM (quantum mechanical–molecular mechanical) simulation method [4]. In this method, the solute molecule is treated quantum mechanically, whereas surrounding solvent molecules are approximated by molecular mechanical potentials. This idea is also used in biological systems by regarding a part of the system, e.g., the activation site region of an enzyme, as a quantum ‘‘solute,’’ which is embedded in the rest of the molecule, which is represented by molecular mechanics. The actual procedure used in this method is very simple: The total energy of the liquid system (or part of a protein) at an instantaneous configuration, generated by a Monte Carlo or molecular dynamics procedure, is evaluated, and the modified Schro¨dinger equations are solved repeatedly until sufficient sampling is accumulated. Since millions of electronic structure calculations are needed for sufficient sampling, the ab initio MO method is usually too slow to be practical in the simulation of chemical or biological systems in solution. Hence a semiempirical theory for electronic structure has been used in these types of simulations.

The QM/MM methods have their own disadvantages, the obvious one being the computational load added to the already complex calculation of the electronic structure.

C. Reference Interaction Site Model

The integral equation method is free of the disadvantages of the continuum model and simulation techniques mentioned in the foregoing, and it gives a microscopic picture of the solvent effect within a reasonable computational time. Since details of the RISM-SCF/ MCSCF method are discussed in the following section we here briefly sketch the reference interaction site model (RISM) theory.

The statistical mechanics of liquids and liquid mixtures has its own long history, but the major breakthrough toward the theory in chemistry was made by Chandler and Andersen in 1971 with the reference interaction site model (RISM) [5]. This theory can be regarded as a natural extension of the Ornstein–Zernike (OZ) equation for simple atomic liquids to a mixture of atoms with chemical bonds represented by intramolecular correlation functions. Introducing this correlation function enables us to take into account the geometry of molecules. However, it cannot handle electrostatics in its original form, though the charge distribution in a molecule plays an essential role in determining the chemical specificity of the molecular system. The next important development in the theory was made in 1981 with the extended RISM theory [6–8]. The extended RISM theory takes into account not only the geometry but also the charge distribution of a molecule, which completes the chemical characterization of a species for the statistical mechanics of a molecular liquid.

A general expression of the RISM equation for a system consisting of several molecular species can be written as

ρhρ ωc ω ω c ρhρ

(1)

where the asterisk indicates convolution integrals and matrix products, and ρ denotes a diagonal matrix consisting of the density of molecular species. c and h are the direct and total correlation matrices with matrix elements of cαβ(r) and hαβ(r), respectively. These functions, cαβ and hαβ, represent the intermolecular correlation functions between the sites (atoms or atomic groups) α and β. ω is the intramolecular correlation function, which

420	Hirata et al.

embodies information about a molecular geometry. In the case of rigid molecules, ω is expressed by

ωαβ(r) ραδαβδ(r) (1 δαβ)	1	δ(r Lαβ)	(2)
	4πLαβ2

where Lαβ indicates the rigid constraint (bond length) representing the chemical bond between sites α and β.

The general equation can be further reduced to the case of infinite dilution limit, a binary mixture, ionic solutions, and so on. These equations are supplemented by closure relations such as the Percus–Yevick (PY) and hypernetted chain (HNC) approximations.

cαβ(r) exp[ βuαβ(r)]{1 hαβ(r) cαβ(r)} 1 (PY)	(3)
and
cαβ(r) exp[ βuαβ(r)hαβ(r) cαβ(r)] {hαβ(r) cαβ(r)} 1	(HNC) (4)

where uαβ(r) is the intermolecular interaction potential and β 1/kBT. All the thermodynamic functions can be calculated from the correlation functions. The excess chemical potential or the solvation energy of a molecule has particular importance and is calculated from the correlation functions [9].

ρ	∫dr cαs(r)	1	2	1		hαs(r)cαs(r)
∆µ β			hαs(r)				(5)
		2		2
	αs

Essentially, the RISM and extended RISM theories can provide information equivalent to that obtained from simulation techniques, namely, thermodynamic properties, microscopic liquid structure, and so on. But it is noteworthy that the computational cost is dramatically reduced by this analytical treatment, which can be combined with the computationally expensive ab initio MO theory. Another aspect of such treatment is the transparent logic that enables phenomena to be understood in terms of statistical mechanics. Many applications have been based on the RISM and extended RISM theories [10,11].

II. OUTLINE OF THE RISM-SCF/MCSCF METHOD

We recently proposed a new method referred to as RISM-SCF/MCSCF based on the ab initio electronic structure theory and the integral equation theory of molecular liquids (RISM). Ten-no et al. [12,13] proposed the original RISM-SCF method in 1993. The basic idea of the method is to replace the reaction field in the continuum models with a microscopic expression in terms of the site–site radial distribution functions between solute and solvent, which can be calculated from the RISM theory. Exploiting the microscopic reaction field, the Fock operator of a molecule in solution can be expressed by

Fsolvi Fi fi bλVλ	(6)
λ solute

where the first term on the right-hand side is the operator for an isolated molecule and the second term represents the solvation effect. bλ is a population operator of solute atoms, and Vλ represents the electrostatic potential of the reaction field at solute atom λ produced

RISM-SCF/MCSCF for Processes in Solutions

421

by solvent molecules. The potential Vλ takes a microscopic expression in terms of the site–site radial distribution functions between solute and solvent:

Vλ solute ρ	qα	gλα(r) dr	(7)
	r
α solvent

where qα is the partial charge on site α, ρ is the bulk density of the solvent, and gλα is the site–site radial distribution function (RDF) or pair correlation function (PCF).

In the RISM-SCF theory, the statistical solvent distribution around the solute is determined by the electronic structure of the solute, whereas the electronic structure of the solute is influenced by the surrounding solvent distribution. Therefore, the ab initio MO calculation and the RISM equation must be solved in a self-consistent manner. It is noted that ‘‘SCF’’ (self-consistent field) applies not only to the electronic structure calculation but to the whole system, e.g., a self-consistent treatment of electronic structure and solvent distribution. The MO part of the method can be readily extended to the more sophisticated levels beyond Hartree–Fock (HF), such as configuration interaction (CI) and coupled cluster (CC).

The solvated Fock operator can be naturally derived from the variational principles [14] defining the Helmholtz free energy of the system ( ) by

ˆ	ˆ	(8)
Ψ|H ∆µ|Ψ Enuc
Here, Enuc is the nuclear repulsion energy and
Ψ|∆µ|Ψ ∆µ		(9)
ˆ

is the modified version of the solvation free energy originally defined by Singer and Chandler [9] [Eq. (5)]. This is a functional of the total correlation function hαs(r), the direct correlation function cαs(r), and the solute wave function |Ψ. Energy of the solute molecule Esolute is defined as follows:

Esolute Ψ|H|Ψ Enuc	(10)
ˆ

Note that this is also a functional of hαs(r), cαs(r), and |Ψ. Imposing constraints concerning the orthonormality of the configuration state function (C) and one-particle orbitals (φi) on the equation, one can derive the Fock operator from based on the variational principle:

	δ( [c, h, t, v, C] [constraints to orthonormality]) 0							(11)
	and
	Fsolvij Fij γij bλ	∂			ρ	∫dr exp[ βuαs(r) hαs(r) cαs(r)]		(12)
		∂qλ
	λ solute				β	αs

If classical Coulombic interactions are assumed among point charges for electrostatic interactions between solute and solvent, and the term for the CI coefficients (C) is omitted, the solvated Fock operator is reduced to Eq. (6). The significance of this definition of the Fock operator from a variational principle is that it enables us to express the analytical first derivative of the free energy with respect to the nuclear coordinate of the solute molecule Rα,