Hideaki Takahashi

Division of Chemical Engineering, Department of Materials Engineering Science, Graduate School of Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, Japan

TABLE OF CONTENTS

1. Abstract

2. Introduction

3. Quantum chemical approach to free energy computation

3.1. Real-space grids method

3.2. QM/MM method combined with a theory of solutions

4. Applications to biological systems

4.1. Electron transfer in biological systems

4.2. Reduction free energy of cofactor in water

4.3. Reduction free energy of cofactor embedded in apoprotein

5. Conclusions and outlook

6. Acknowledgement

7. References

1. ABSTRACT

We recently developed a novel computational methodology, referred to as QM/MM-ER, to compute free energy change associated with a chemical reaction in a condensed phase by combining the quantum mechanical / molecular mechanical (QM/MM) method with a theory of solutions. In the present review article we illustrate an outline of the QM/MM-ER method. We also present an extension of the QM/MM-ER method to compute reduction free energy of cofactor FAD (Flavin Adenine Dinucleotide) in water as well as in apoprotein by introducing a novel approach. The key of the approach is that only the excess charge involved in the reduction process is identified as a solute. The adequacy of the method is examined for the reaction in aqueous solution by comparing the result with that obtained by a conventional approach. The reduction free energy of cofactor FAD embedded in a cholesterol oxidase (PDB id = 1B4V) is computed as 163.1 kcal/mol by the QM/MM-ER approach, while it is obtained as 80.1 kcal/mol for the cofactor in water solution.

2. INTRODUCTION

The quantum chemical calculation (1) has been extensively applied to various systems to study molecular structures, properties, and reactivities. So far the efficiency of the method has been well established. The major origin of the success of the quantum chemistry lies in the abilities that it provides us with reliable data with substantial accuracies and also in the fact that it enables us to analyze the behavior of the electrons involved in the system. Indeed, the source of physical properties of interests as well as the underlying mechanisms for chemical reactions can be successfully interpreted in terms of the molecular orbitals (MO) which provides a useful picture based on the approximation that individual electrons are moving independently. The Hartree-Fock (HF) approach in the MO theory is constructed on the basis of the mean-field approximation and neglects the dynamic correlation between electrons. Within more sophisticated approaches such as configuration interaction (CI) one can take into account the electron correlations by involving the electronically excited configurations generated from the HF reference state. Such a procedure allows a systematic improvement on the HF wavefunction. However, the computational cost increases drastically as the number of electrons involved in the correlation increases. Hence, the applications of the MO theories had been limited only to a gas phase reaction with a small number of electrons, although most of the chemical reactions take place in a condensed phase such as solution. During the past several decades, there were great strides in the method of the electronic structure calculations. The density functional theory (DFT) (2), specifically the Kohn-Sham DFT (3, 4), enables ones to perform the first-principles calculations for relatively larger systems. It achieved a great success in describing the total energy as well as the electron density of the system in the ground state by projecting the non-local exchange and correlation potential onto a local potential as a functional of the electron density. Indeed it enables ones to compute energy of the system in comparable accuracies with the time-consuming MO theories within much less computational costs. There have also been a lot of theoretical developments to refine the exchange-correlation functionals which dominate the quality of the DFT calculation. Due mainly to the success in DFT, the quantum chemical approach is now being applied to various biological systems.

A biological system including the surrounding water solvent is composed of many particles, hence, we have to make a device to treat quantum chemically such a large system. The Car-Parrinello (CP) approach (5) opened the way to the first-principles molecular dynamics simulations for reacting condensed systems. They circumvented the explicit diagonalization of the electronic Hamiltonian at each time step of the molecular dynamics simulation by propagating the Kohn-Sham wavefunctions in the same time scale with the nuclear motion, for which a fictitious mass for the electron is introduced. During the Car-Parrinello simulation, the electronic wavefunction of the system is being loosely bound to the solution given by the Born-Oppenheimer approximation. The CP method is efficient, however, it necessitates rather small time step in the numerical time propagation as compared with the classical molecular dynamics simulations, and it still needs a massively parallel computer to execute a simulation that involves a few hundreds of atoms. It is often the case that only a small part of the whole system participates in the chemical event and the remaining part of the system serves as a static environment. The hybrid quantum mechanical/molecular mechanical (QM/MM) method (6, 7) provides an appropriate framework to handle such a system by dividing the whole system into two subsystems. The chemically active site of the system is described by quantum chemistry, while the electronically static environment is represented by a classical force field. The effect of the electrostatic field formed by a set of point charges in the MM subsystem on the reaction is taken into consideration by putting the electrostatic field in the Hamiltonian of the QM subsystem. So far the QM/MM method has been widely applied to various systems such as solutions or biological systems, and its efficiency has been well established. Recently, we developed a novel QM/MM approach (8-12) that utilizes real-space grids (13, 14) to express the electronic wavefunctions in the QM subsystems described by Kohn-Sham DFT. The advantage of the use of the real-space grids approach is that the parallel implementation is amenable due to the fact that the effective Hamiltonian in the real-space representation is almost diagonal within the framework of the Kohn-Sham DFT. In this article, we review the QM/MM method with real-space grids and also discuss the parallel efficiency with the method.

The quantity that plays a central role in determining the reaction path in the many-particle system with the thermal fluctuation is, of course, the free energy change associated with the reaction. Hence, the computation of the free energy is definitely an issue of crucial importance in the study of the chemical reaction in the condensed system. However, the evaluation of the free energy with the method of molecular simulation is known to be a very demanding task. The free energy perturbation (FEP) or thermodynamic integration (TI) method is based on the Kirkwood charging formula (15-17) and it introduces a set of intermediate points on an arbitrary path connecting the initial and the final states of a chemical process of interest. Then, the free energy changes between the adjacent points are accumulated along the path to evaluate the free energy change associated with the reaction. The FEP or TI is numerically exact, however, it should be noted that a sufficient number of intermediate points, which are physically of no importance, must be prepared to avoid being ill-sampled for the molecular configurations. Also it is required for the convergence in the free energy to sample a large number of molecular configurations at each intermediate point. Therefore, the free energy calculation by the FEP approach along with the quantum chemical calculation becomes almost intractable due to the twofold time-consuming procedures. The quest is, thus, required for an efficient quantum chemical approach with which one can compute free energy change under the influence of the environment.

The polarizable continuum model (PCM) approximates the solvent as a continuum with dielectric constant (18). In the field of quantum chemistry the PCM approach is extensively utilized due mainly to its practical efficiency. The deficiency of the approach is that it completely neglects the structure of the constituent molecules of a solvent and cannot represent the short range interactions such as hydrogen bonds. In the context of the density functional theory of solutions, the solvation free energy of a solute molecule can be described in terms of the distribution function of the solvent around the solute. In general, the position of the solvent molecule with respect to the solute is specified by a multi-dimensional coordinate. In the practical implementation, however, the method of reference interaction site model (RISM) is often used (19, 20), where the solvent distribution is approximated by a set of site-site radial distribution functions. An approximate set of integral equations is, then, introduced to obtain the site-site radial distribution functions by adopting closure relations such as PY or HNC approximations (15). A remarkable advance, referred to as RISM-SCF (21), was made by combining the quantum chemical SCF (self-consistent field) calculation with the RISM integral equations. An important aspect of the RISM-SCF approach is that it soundly incorporates the effect of the solvent as a structured environment into the quantum chemical calculation. This is definitely more advantageous than the PCM approach. The RISM-SCF method, however, has a drawback that the electron density of the solute molecule is to be reduced to a set of point charges placed on the sites in the procedure to construct radial distribution functions. It indicates that the cloud-like electron density of the solute, that is an inherent nature of the quantum mechanical object, is inevitably spoiled in the description of the solvation free energy. Indeed, our recent study revealed that the neglect of the volume of the electron density leads to erroneous results in the computation of the free energy difference due to the overestimation of the solute-solvent interactions especially at the anionic site. It is also well known that RISM-SCF cannot be successfully applied to the low density region of solution.

In this article, we review our new quantum chemical method (22) recently developed to compute free energy change for chemical reactions in condensed system. Our strategy to overcome the difficulty in the free energy calculation is to employ the hybrid quantum mechanical / molecular mechanical (QM/MM) approach with the real-space grids in combination with a novel theory of solution, the theory of energy representation (QM/MM-ER). In the method of the energy representation (23-25), the solvent distribution is constructed with respect to the solute-solvent interaction potential and it serves as a fundamental variable to describe the free energy instead of the spatial distribution functions. An emphasis must be placed on the fact that the concept of the interaction site is no longer needed in the construction of the energy distribution functions and, therefore, the combination with the quantum chemical object is straightforward. The diffuseness of the electron density of the solute and its fluctuation can be naturally incorporated into the free energy calculations without introducing special devices for it. In the previous works (22, 26-30), we applied the method to various chemical events in aqueous solutions to examine the efficiency and the reliability of the method and found that it reproduces the solvation free energies in excellent agreement with the experimental results. The extension of the method is now in progress for the chemical reactions in biological systems, such as ligand binding at active site, redox reaction in the cofactor, proton binding, and the chemical bond rearrangement. In the following, we make a review for the novel methodology and its application to the computation of the reduction free energy of a cofactor which is embedded both in water and in an apoprotein.

3. QUANTUM CHEMICAL APPROACH TO FREE ENERGY COMPUTATION

In our recent development, referred to as QM/MM-ER (22), we combined the QM/MM method with the theory of energy representation to compute free energy change in condensed systems. We described the electronic structure of the QM subsystem by the Kohn-Sham DFT, where the one-electron wavefunctions are represented by real-space grids to achieve high parallel efficiency. In Subsection 3.1 at first we make a brief review for the real-space grids approach placing a focus on the parallel implementation and its efficiency. In Subsection 3.2 we illustrate the formulation to combine the method of QM/MM with the theory of the energy representation.

3.1. Real-space grids method

In the field of quantum chemistry, of which major target is atoms or molecules, the LCAO (linear combination of atomic orbitals) method has been extensively utilized to expand the one-electron wavefunction in the numerical implementation of the molecular orbitals (MO) and the density functional theory (DFT). It is known that the Slater-type orbital is more appropriate than the Gaussian orbital as a component atomic orbital to reproduce the realistic behaviors of molecular orbitals especially near the atomic core regions. However, the Gaussian-type orbitals are commonly used instead of Slater-type ones because of the numerical convenience. In practice, a basis function used in the LCAO approach is constructed from a linear combination of primitive Gaussians with fixed contraction coefficients so that it mimics the behavior of a Slater-type orbital. As a result, a huge number of two electron integrals with respect to primitive basis functions are to be computed and stored in preparation for the self-consistent field (SCF) calculations. The adequacy of the LCAO approach is well established by a number of numerical applications. A critical drawback in the LCAO approach is, however, that the effective Hamiltonian represented by atomic orbitals has non-zero elements in its off-diagonal part, which gives rise to a large amount of data communications among the processors in the parallel implementations. In the study for the electronic structures of infinite systems such as crystals, on the other hand, plane wave basis set is often employed to make use of its periodicity. The matrix for the kinetic energy operator in the effective Hamiltonian becomes diagonal in the plane wave basis representation, which is obviously advantageous for the parallel computation. However, other operators in the Hamiltonian are not diagonal in the momentum space representation. Moreover, the approximate exchange-correlation potential which appears in the Kohn-Sham equation can be evaluated only after the total electron density is obtained in the real space. It is, therefore, required to transform one-electron wavefunctions represented in the momentum space to those in the real-space, where the fast Fourier transformation (FFT) algorithm is often utilized to expedite the transformations. The FFT is also needed for the backward transformation from real to momentum space for each one-electron wavefunction. The computational cost of FFT for an orbital in general scales as Nlog₂N where N is the number of plane waves used in the computation and FFT is the most time-consuming part in the method of the plane wave basis. It should also be noted that FFT is never suitable for the parallel computing especially when one uses a parallel computer with a distributed memory architecture.

In a recent development Chelikowsky proposed a method of constructing a one-electron wavefunction on the real-space grids (31-33). Since most of the operators in the Kohn-Sham equation is local in the real-space representation it is quite natural to express the one-electron orbitals by a set of probability amplitudes defined on the discrete grid points that are uniformly distributed over a real-space cell. In the parallel implementation, we distribute the values for the wavefunctions as well as the operators to the processors by dividing the cell into subspaces. This procedure is schematically illustrated in Fig. 1. The Hamiltonian matrix has non-zero elements only in the vicinity of diagonal part, and hence, the data-communication among the processors can be suppressed minimally, which benefits the high performance in the parallel computation. Moreover, the real-space grids approach has several advantages as follows. At first the local augmentation of the basis can be done straightforwardly by introducing dense grids around atomic cores for instance. Second, the periodic and the non-periodic boundary conditions are available without serious modifications of the code. Third, the overlap matrix of the real-space grids basis becomes an identity matrix and no additional procedure is needed to orthogonalize the basis set. It should be noted, however, that the drawbacks also exist in the real-space approach. The electronic energy of an atom fluctuates for the variation of the relative position of the atom with respect to the grid points. Such a phenomenon is known as the "egg-box problem" (13) and it can be substantially alleviated by the use of the dense grids near the atomic core regions (14). In addition the molecular integral evaluated in the real-space approach is less accurate as compared to the atomic orbital basis. The most difficult problem in the real-space method arises from the evaluation of the operator with serious non-locality such as Hartree-Fock exchange operator. The real-space representation is of great advantage to the local operator, however, it becomes intractable to compute integrals associated with non local operators. Fortunately, most of the operators in Kohn-Sham DFT are local or semi local except for the hybrid type functional which contains the Hartree-Fock exchange. In the following paragraphs we describe the real-space approach for the implementation of the Kohn-Sham DFT. An emphasis will be placed on the efficiency of the parallel computations.

In the framework of the Kohn-Sham DFT, the Schrödinger equations for the one-electron orbitals {_i(r)}, which constitute a non-interacting reference system, are given in the atomic units as

(3.1)

where i denotes the suffix allotted to each orbital, n(r) is the total electron density at position r, and N is the number of occupied orbitals. The terms in the square bracket in Eq. (3.1) indicate in order the kinetic energy operator, Hartree potential, electron-nuclei pseudopotential and the exchange-correlation potential for electrons. In the following we discuss for each term the applicability of the real-space representation to the parallel implementation. The kinetic energy operator is expressed by the finite-difference scheme in the real-space representation. Here, we formulate the matrix elements of the operator with respect to a one-dimensional wavefunction for simplicity. The Taylor expansion of a wavefunction at grid point l is given as

(3.2)

where h is the grid spacing. By making the sum of upper and lower equations of Eq. (3.2), we obtain the first order finite difference expression for the second order derivative , thus,

(3.3)

The higher order finite difference expression can also be formulated by considering the higher order terms in the Taylor expansion of Eq. (3.2). Then, the general form of the kinetic energy operator in 3 dimension is represented as,

(3.4)

where L denotes the order of the expansion and C_k1, C _k2, and C _k3 are the expansion coefficients for which a Table is given by Chelikowsky for L = 1 ~ 6 (32). According to our preliminary calculations it is found that adopting L = 4 is sufficient enough to achieve a desired accuracy. In Eq. (3.4) we can readily recognize that the non-locality of the kinetic energy operator in the finite-difference representation is confined within a small region of space. Eq. (3.4) suggests that the operation on the wavefunction at the grid point necessitates the values of the wavefunction only on the neighboring 8 grid points for each direction in the case of the 4th order finite-difference method (L = 4).

The second term in the square bracket in Eq. (3.1) is termed as the Hartree potential v_H(r) and it represents the classical coulomb potential formed by the electron density n(r). In the periodic system, v_H(r) can be evaluated by transforming the electron density in real space to the reciprocal space by using FFT. In the non-periodic system, on the other hand, v_H(r) can be determined by the direct integration of the classical electron repulsions between the electron densities on all the grid point. However, it is computationally very demanding since it requires two-fold loop that runs over all the pair of grid points in the real-space cell. Instead, we adopt a method solving the Poisson equation for the Hartree potential, thus,

,(3.5)

for which the conjugate gradient (CG) algorithm synchronized with the steepest descent (SD) procedure is employed (14). Eq. (3.5) involves the second-order derivative of the Hartree potential, hence the amount of data to be communicated among the processors is same as that needed for the kinetic energy operator. We also note that multicenter numerical integration scheme (34) proposed by Becke is used to determine the boundary condition for the solution of Eq. (3.5).

The third term in Eq. (3.1) is the pseudopotential (35,36) formed by the core electrons and the nuclei, that is introduced to realize the smooth behavior of the wavefunction. The use of the pseudopotential is essential for the real-space grid approach since the true wavefunction oscillates rapidly near the atomic core regions. In the present calculation we employ the norm-conserving pseudopotential with Kleinman-Bylander separable form (36). It includes the non-local part and the overlap integrals must be evaluated between the Kohn-Sham orbitals and the atomic pseudopotentials. However, it can also be incorporated easily into the parallel computation without serious loss of efficiency. This is due to the fact that the non-local pseudopotential has non-zero values only within a small region of space near the atomic core. To improve the description of the rapid behavior of the non-local part of the pseudopotential, we employ the time-saving double grid method developed by Ono and Hirose (37). The feature of the method is to neglect the explicit computation of the Kohn-Sham orbitals on dense-grid points and they are estimated by Lagrange interpolations of the original coarse-grid points instead. The interpolation can be done successfully by the benefit of the smooth behavior of the wave functions on the pseudopotentials. The contribution of the dense-grid points is implicitly involved in the weight factors defined at coarse-grid points. It is worthy of note that there is no need to recompute the weight factors on coarse grid points during the SCF procedure. The accuracy of the real-space grid approach reinforced by the time-saving double-grid method has been well established by applying the method to various systems. The detailed explanation of the double grid method is presented in Ref. 37. The computational accuracy is almost comparable to that of the sophisticated Gaussian basis sets such as Dunning's correlation consistent basis set with polarization augmented by diffuse functions (aug-cc-pVDZ) (38).

The last term to be examined is the exchange-correlation potential v_xc(r) in Eq. (3.1). A rude estimation of the exchange-correlation energy E_xc is usually given by the local density approximations (LDA) (2) which assumes that E_xc can be expressed as

.(3.6)

 in Eq. (3.6) denotes the exchange-correlation energy per electron in a uniform electron gas with density (r). The exact form of the exchange part E_x as a functional of (r) was given by Dirac in 1930 for the uniform electron gas (39). Of course, Eq. (3.6) is valid only for the system of which electron density gradually changes. Although realistic atoms or molecules are far from this situation in general, it is known that Eq. (3.6) is adequate at least for qualitative discussions. In a sophisticated approach, the LDA is corrected by a functional of the gradient of the electron density, which is referred to as the generalized gradient approximation (GGA). We apply the gradient correction scheme for exchange energy E_x proposed by Becke (40), which is expressed as

,(3.7)

where parameter  is set at 0.0042, is the electron density for spin  and x_is termed as a inhomogeneity parameter which is defined as

.(3.8)

Note that the gradient for the electron density with spin  appears in Eq. (3.8). The above GGA functional or its hybridization with the exact Hartree-Fock exchange (41) is known to be comparable in accuracy to the sophisticated methods in MO theories. The point to be noted here is that the exchange functional in the form of Eq. (3.7) only necessitates the gradient of the electron density at most. It suggests that the information contents to be communicated in the parallel computation are common to those required in the operation for kinetic energy. The situation is almost the same with the correlation energy E_c formulated by Lee, Yang, and Parr (LYP) (42). Thus, the parallelization of the exchange-correlation potential is straightforward even at the GGA level.

Thus, it has been found that the operators in the effective Hamiltonian for the Kohn-Sham DFT are almost local or semi local in the real-space representation. Although it includes the non-local operators, its non-locality is limited only in the small region of space. Hence, it is expected that the real-space grids approach can afford a high performance in the parallel computing. We performed a benchmark test calculation for a system containing 98 valence electrons by using our original optimized code. The data communication is commanded by the standard parallel interface MPI. The parallel efficiency of the real-space grid algorithm is obtained as 80.5 % on a computer cluster with 8 CPUs (Pentium 4/3.8 GHz) connected by ordinary communication facilities (PCI/Giga bit Ethernet). This demonstrates the efficiency of the real-space grids approach in combination with the Kohn-Sham DFT.

3.2. QM/MM method combined with a theory of solutions

In most cases, only a small part of a system takes part in the chemical event of interest in the condensed phase such as a solution or a biological system. The hybrid quantum mechanical/molecular mechanical (QM/MM) approach provides a reliable computational framework to describe such a situation, where the chemically active solute is described by quantum chemistry and the remaining static environment is represented by a classical empirical force field (16,17). Within the framework of the QM/MM method, the total energy E_T of the system is decomposed into three terms, thus,

(3.9)

where E_QM is the electronic energy of the QM subsystem including the nuclear repulsion energy, E_MM is the energy of the MM subsystem represented by classical force field, and the term E_QM/MM denotes the interaction between the QM and the MM subsystems. The electronic wavefunction  with an eigenvalue E under an instantaneous solvent configuration can be defined by the following Schrödinger equation,

(3.10)

where H₀ is the electronic Hamiltonian of the isolated QM subsystem including the nuclear repulsion energies while V_ELS is the electrostatic field formed by the point charges placed on the interaction sites of the MM molecules. Then, we can define the energy E_QM in Eq. (3.9) as

.(3.11)

Here, we introduce the energy termed as distortion energy E_dist of the QM subsystem, thus,

(3.12)

where E₀ and _ are the energy and the wavefunction of the isolated solute, respectively. The interaction energy E_QM/MM in Eq. (3.9) is expressed by

(3.13)

where E_ELS is the Coulombic interaction between the wavefunction and the MM molecules and given by

.(3.14)

In Eq. (3.13) E_ns is the Coulombic interaction between the nuclei in the QM subsystem and the point charges on MM molecules, and E_vdW is the van der Waals interaction between the QM and the MM subsystems. The last term in Eq. (3.9) is the energy of the MM molecules described in terms of the two-body interactions defined by a classical force field. In our implementation of the QM/MM method, we adopt the Kohn-Sham DFT utilizing the real-space grids described in the previous subsection.

The solvation free energy is of particular importance for the study of chemical reaction in solution since the reaction free energy in solution can be described in terms of the solvation free energies of the initial and the final states of the chemical event and the free energy change in the gas phase. Therefore, the most important quantity to know the reaction path in the condensed phase is obviously the solvation free energies of the relevant molecules. The solvation free energy  of a QM solute can be expressed as

(3.15)

where  is the inverse of the product of the Boltzmann constant k_B and temperature T, and n is the shorthand notation of the instantaneous electron density n(r) of the QM solute under a given solvent configuration denoted collectively as X. Recently we developed a methodology to compute solvation free energy  in the form of Eq. (3.15) by combining the QM/MM procedure with the theory of energy representation (QM/MM-ER). Within the theory of energy representation (23-25), the solvation free energy is exactly formulated in terms of the distribution functions of the solute-solvent interaction potential. A characteristic feature of the method of the energy representation is that it no longer needs the concept of the interaction sites. Hence, the diffuse nature of the electron density in the quantum chemical object can be naturally incorporated into the free energy calculation. However, since the standard version of the method of energy representation assumes that the solute-solvent interaction is pairwise additive, a treatment must be made in combining the method with the QM/MM approach that involves many-body interactions. To overcome this difficulty we resort to a two-step approach, in which the free energy contribution due to the two-body interaction is evaluated in the standard framework of the energy representation and the remaining minor contribution is computed separately. The key of this approach is to construct an intermediate state of the QM solute with the electron density fixed at its average distribution in solution. We review the outline of the methodology in the following. The free energy  can be decomposed into three terms, thus,

(3.16)

where is the average distortion energy E_dist defined by Eq. (3.12), is the solvation free energy of a QM solute with the average electron density in solution, and  is the free energy due to the electron density fluctuation around its average distribution . and are simply statistical averages and, respectively, given by

,(3.17)

.(3.18)

Note that Eqs. (3.17) and (3.18) are easily obtained through an ordinary QM/MM simulation simultaneously. The solvation free energy can be regarded as a contribution due to the two-body interaction between the QM solute and the solvent and is expressed by

,(3.19)

while the remaining minor contribution  due to the many body effect of the QM object can be given by

.(3.20)

It can be readily understood that Eq. (3.15) is reproduced by substituting Eqs. (3.19), and (3.20) to Eq. (3.16). It should also be noted that Eq. (3.16) holds independently of and under the definitions of Eqs. (3.19) and (3.20). We take the choices for and as defined by Eqs. (3.17) and (3.18) for the conceptual and the numerical convenience. In the computation of the free energy , the solute-solvent interaction is pairwise additive since the electron density of the solute is fixed at its average distribution and no longer corresponds to the instantaneous solvent configurations. Therefore, the theory of energy representation can be directly applied to compute the free energy . The remaining part of the free energy due to the electron density fluctuation can also be evaluated separately within the framework of the energy representation. Next, we make a brief review for the procedure to obtain along with the explanation for the theory of energy representation.

In the theory of energy representation, the energy distribution functions for the solvent, instead of the spatial distribution functions, serve as fundamental variables to determine the free energy change of interest. At first we define the instantaneous energy distribution function as

(3.21)

where x_i is the full coordinate specifying the position and the orientation of the ith solvent and v is a potential function and termed as the defining potential introduced to define the energy coordinate of the solvent molecules. Usually v is taken as the solute-solvent interaction potential under consideration. Then, the energy distribution function _u() for the solvent molecules under the condition that the solute-solvent interaction is u is defined as the ensemble average of the instantaneous energy distribution, thus,

(3.22)

where U stands for the solvent-solvent potential energy and X specifies the configuration of the solvent. When we take u = 0 in Eq. (3.22) the system is termed as a reference (or a pure solvent) system and the resulting energy distribution function is denoted as ₀(). On the other hand, when u is taken as the solute-solvent interaction potential v, the system is identified as a solution system and the distribution function is simply written as () in particular. The solvation free energy defined by Eq. (3.19) can be exactly formulated in terms of the energy distribution functions, thus,

.(3.23)

In Eq. (3.23), the energy distribution functions ₀() and () are constructed with respect to the solute-solvent interaction energy between the QM solute with electron density and the solvent molecules for the purpose to obtain . In Eq. (3.23) () is the indirect part of the solute-solvent potential of mean force and  is the coupling parameter associated with the gradual insertion of the solute into the solvent. By virtue of the one-to-one correspondence between the set of interaction potentials and the resulting energy distribution functions, () can be formulated by a functional of the energy distribution function. In the practical implementation, integration of () with respect to  is evaluated by adopting a hybrid functional of PY and HNC approximations (15). We refer the readers to the previous works for a explicit form of the functional (22,24). The free energy change  of Eq. (3.20) that originates from the many-body interaction involved in the QM solute can also be computed in terms of the energy distribution functions constructed over the energy coordinate of . In Fig. 2 a flow chart is presented to show the procedure to compute .

(Step 0) Classical molecular dynamics simulation is performed to relax the configuration of the whole system.

(Step 1)An ordinary QM/MM simulation is carried out to obtain the average distortion energy and which are defined by Eq. (3.17) and (3.18), respectively.

(Step 2)The energy distribution functions ₀(), ₀() and () for pure solvent and solution systems are computed for the QM solute with the electron density fixed at the distribution of .

(Step 3)For the computation of the free energy contribution  due to the electron density fluctuation, additional QM/MM simulations have to be performed. Then, the energy distribution functions , and () for pure solvent and solution systems are computed with respect to the energy coordinate  defined above.

For the definition of the two dimensional matrices ₀() and in the flow chart we refer the readers to Ref. 22.

There are notable advantages in the method of QM/MM-ER in the computation of the free energy changes in condensed phases. By virtue of the theory of energy representation, which treats the solvent molecules with the same energy coordinate as an identical object, the convergence of the distribution function on the energy coordinate is much faster as compared with that of the spatial distribution function. This benefits substantially the computation of the free energy change based on the quantum chemistry since the quantum chemical calculation requires a large computational costs in general. Since the energy distribution function is obtained by the projection of the spatial distribution function onto the energy coordinate, the information content is substantially reduced indeed. It is worthy of note, however, that the free energy change can be formulated without approximation in principle. In addition, the diffuseness of the electron density as well as its fluctuation in solution can be naturally incorporated into the free energy calculation. Our previous work suggests that the neglect of the diffuse nature of the electron density may lead to erroneous results in some cases. The accuracy and the efficiency of QM/MM-ER approach was well examined in the previous works. For instance in Ref. 8 and 22 the solvation of a water molecule in water solution was investigated by the real-space grid QM/MM approach combined with the theory of energy representation and it was demonstrated that the solvation free energies as well as the solvation structures were accurately produced. In Ref. 26 the pKw values in the supercritical states of water, which are difficult to be determined experimentally, were successfully computed. In Ref. 27 the free energy change was obtained with substantial accuracy for the isomerization reaction of glycine in aqueous solution and a quantitative discussion was first given about the importance of the inclusion of the diffuseness of the electron density at an anionic site in the description of the solute-solvent interactions. In Ref. 28 it was revealed that the QM/MM-ER approach is adequate enough to predict even the slight free energy difference associated with the anti/syn conformation change for acetic acid in aqueous solution. Further, in Ref. 29 the relative acidity of acetic acid with respect to water was determined in good agreement with the experimental value. Finally, in Ref. 30 a novel reaction path was suggested by performing the QM/MM-ER simulations, where a chemical reaction in hot water without catalyst is substantially promoted by the proton transfers along the hydrogen bonds of water molecules. We refer the readers to these articles for more details.

4. APPLICATIONS TO BIOLOGICAL SYSTEMS

So far the QM/MM-ER method has been applied to various chemical reactions in solution and its efficiency and accuracy have been established. However, we may encounter some difficulty in extending the method to a chemical event taking place in a part of a protein immersed in water. Since the solvation free energy of a protein amounts to a large negative value, it is numerically problematic to evaluate the reaction free energy in terms of the solvation free energies of the proteins in the initial and the final states of the reaction. In this Section we propose a methodology to compute free energy change associated with an electron transfer (i.e. oxidation and reduction) reaction in solution or protein. To circumvent the problem noted above we take a novel approach where only the excess charge involved in the redox reaction is regarded as a solute. In Subsection 4.1, we introduce a redox system (cofactor FAD) chosen as a test calculation. In Subsection 4.2, the outline of the methodology is illustrated for the computation of the reduction free energy of a cofactor immersed in water within the framework of QM/MM-ER. The validity of the novel approach is also examined by performing a conventional method. In Subsection 4.3 the method is extended to the free energy calculation for the reduction of the cofactor embedded in the apoenzyme.

4.1. Electron transfer in biological systems

The electron transfer, i.e. oxidation and reduction, is a ubiquitous elementary reaction among the variety of chemical processes in biological systems. Indeed, it plays a central role for instance in storing the free energy used for the ATP synthesis by pumping protons against the gradient of the proton concentration (43). Hence, it is an intriguing problem how the biological environment regulates the electron translocations. Since the fundamental quantity that governs the reaction path in a condensed phase is the free energy change associated with the process, it is necessitated to establish a methodology to compute free energy for the electron transfer under the influence of a biological environment. In this subsection, we present an illustration for a novel method to compute free energy change for the reduction of a protein immersed in water within the framework of the QM/MM-ER approach described above. The point of the method is to treat the excess electron to be attached on the molecule as a solute in the construction of the distribution functions with respect to the solute-solvent interaction energy. Then, the remaining molecules are considered as solvent.

The system we choose for the test calculation is the cholesterol oxidase (PDB id: 1B4V) in which a redox cofactor FAD (flavin adenine dinucleotide) is noncovalently attached to the protein through hydrogen bonds. FAD acts as an active site of the protein and undergoes one- or two- electron reduction. In the present calculation, we focus our attention on the one-electron reduction and compute the free energy change for it. We also perform the QM/MM-ER simulation for the system where only the active site of FAD, isoalloxazine ring, is immersed in water (44). A comparison will be made to clarify the role of the environmental protein in promoting the electron transfers. In the previous work Q. Cui carried out the QM/MM free energy perturbation calculations for the reduction free energy of the same system (45). A notable feature of their approach is that they employed a coupling scheme termed the dual-topology-single-coordinates scheme, in which the initial (oxidized) and the final (reduced) states of the solute are forced to adopt the same molecular geometry during the perturbation simulation. They circumvented the problem associated with the treatment of the QM object with fractional electrons appearing in the intermediate points between the oxidized and reduced states by introducing separate copies of the two chemical states in the total energy potential. Their approach is numerically rigorous, however, the computational effort is twice as large as the conventional QM/MM free energy perturbation scheme. In the QM/MM simulations, they resorted to the self-consistent-charge density-functional-tight-binding (SCC-DFTB) method (46) in which approximations are made to the two-electron integrals to expedite the quantum chemical calculations. In contrast to their approach, our present method requires the QM/MM simulations only for the initial and the final states of the reduction process and, hence, the computational cost is much more modest.

4.2. Reduction free energy of FAD in water (44)

7,8-dimethyl isoalloxazine ring is the active site of cofactor FAD (Fig. 3(a)). In the present application two methyl groups in the isoalloxazine ring are replaced by two hydrogen atoms as shown in Fig. 3(b). In the reduction process the excess charge is considered to occupy the  orbital of the isoalloxazine ring. Obviously, investigation of the redox property of the isoalloxazine ring is the key step to understand the mechanism responsible for the electron transfer occurring between the substrate and the cofactor. At first we computed the reduction free energy for the isoalloxazine ring immersed in water as a pilot calculation that utilizes the QM/MM-ER method. In this calculation, we test a new approach where the excess charge is regarded as a solute and the rest of the system is treated as solvent in the implementation of QM/MM-ER. In the following, we present an outline of the methodology with some formulations.

In applying the QM/MM approach to the reduction process in water, the isoalloxazine ring is described quantum chemically, while the solvent water molecule is represented by a classical model. Then, the reduction free energy can be expressed in terms of the energy components in Eq. (3.9), thus,

(4.1)

where the notations of N+1 and N specify the systems in the reduced and the oxidized states, respectively. In Eq. (4.1) the electron density of the isoalloxazine ring is denoted by n and the instantaneous molecular configuration of the MM subsystem is collectively expressed as X. To adopt the method of the energy representation for the computation of , only the excess charge is treated as a solute and the isoalloxazine ring and water molecules are regarded as a mixed solvent. Here, we introduce the energy coordinates for both two kinds of solvent with respect to the solute, i.e. excess charge. The energy coordinate for the QM object (isoalloxazine ring) can be introduced by defining the interaction potential between the excess charge and the isoalloxazine ring as

(4.2)

where stands for the electronic wavefunction and H₀ represents the electronic Hamiltonian including the nuclear repulsion energies in the absence of the electrostatic potential of the surrounding water solvent. Note that is obtained as the eigenfunction for the following Schrödinger equation,

.(4.3)

In Eq. (4.3) V_pc is the electrostatic potential formed by the point charges placed on the interaction sites of the water molecules of which full coordinates are collectively represented by X. The equality also holds for the equation obtained by replacing the notations N+1 in Eq. (4.3) by N. The interaction potential v_QM in Eq. (4.2) is interpreted as the energy difference between the electronic states for the reduced and the oxidized solutes excluding the stabilization energies due to the hydrations. Next we introduce the energy coordinate for a water molecule. The interaction potential between the excess charge and ith water molecule is given by

(4.4)

where q_i,m denotes the point charge in the unit of elementary charge of mth site placed on ith water molecule and r_i,m is the position vector of the site. n_excess(r) in Eq. (4.4) is the electron density of the excess charge attached on the isoalloxazine ring and is defined as the difference between the electron densities derived from the wavefunctions and . The energy distribution functions defined by Eq. (3.22) are constructed for both the interaction energies v_QM and defined by Eqs. (4.2) and (4.4), respectively.

Here we describe in detail the reference and the solution systems for which the energy distribution functions are constructed. In the novel approach where only the excess charge is regarded as a solute, the reference system is identified as the oxidized isoalloxazine ring and the surrounding solvent water molecules, while the solution system is composed of the reduced ring and the water molecules. In the reference system, solvent configurations are sampled under the condition that the excess charge is absent in the computation of the interaction between the isoalloxazine ring and water molecules. On the other hand, the interaction between the water molecules and the excess charge is taken into consideration in generating the configurations of the water molecules in the solution system. The individual energy distribution functions _QM() and _MM() for the QM and MM molecules are separately being constructed during a QM/MM simulation. The total energy distribution function () for the mixed solvent system can be given as the sum of these distributions, thus, () = _QM() + _MM(). The same is true for the energy distribution ₀() in the reference system. Then, the reduction free energy can be evaluated by substituting these distribution functions into Eq. (3.23). Of course, the individual free energy contributions from the QM or the MM solvents can also be obtained formally by using the corresponding energy distribution functions for these solvent molecules. It should also be noted that the correlation matrix used in the functional to evaluate the integral of (; ) in Eq. (3.23) is not necessitated since the excess charge does not posses the effective exclusion volume. This leads to the fast convergence in the free energy since the two-dimensional correlation matrix is much slower in convergence than the one-dimensional energy distribution function.

Computational set up for the QM/MM-ER simulation is as follows. The molecular geometries for the QM molecules, i.e. the reduced and the oxidized isoalloxazine rings, are optimized by the Gaussian 03 package (47) with the hybrid UB3LYP functional (41,42) and the aug-cc-pVDZ basis set (38). The QM molecule is placed at the center of a water solvent confined within a spherical cavity with a van der Waals wall. The water droplet consists of 676 TIP3P water molecules (48). In the QM/MM simulations the electronic structure of the QM subsystem is determined by the Kohn-Sham DFT that utilizes the real-space grids. The grid spacing is set at 0.314 a.u. The exchange and correlation energy of electrons are approximately evaluated by the BLYP functional. The fourth-order-finite difference scheme (N = 4) in Eq. (3.4) is adopted for the kinetic energy operator of electrons. Interactions between the valence electrons and the nuclei are represented by the norm-conserving pseudopotentials in Kleinman-Bylander form (36). Molecular configurations of the reduced and the oxidized isoalloxazine rings are fixed during the QM/MM simulations. The Lennard-Jones (LJ) parameters in CHARMM 22 force field (49) are assigned to the atoms in the QM molecule to compute the LJ interactions between the QM and the MM subsystems. The molecular configurations of the water solvent are sampled by performing the molecular dynamics simulations where the Newtonian equations of motions for the water molecules are solved numerically by using the Verlet algorithm (16). To obtain ensemble averages for the energy distribution functions given by Eq. (3.22), 50-ps QM/MM simulations are performed for the solution and the reference systems. The energy distribution functions for water solvent are constructed for the water molecules within a sphere with radius of 15 Å centered at the center of mass of the isoalloxazine ring. The thermodynamic condition of water solution is set at T = 300 K and  = 1.0 g/cm³. In the next paragraphs we present the results for the QM/MM-ER simulations performed for the water solution.

The one-electron reduction energy of isoalloxazine ring at the gas phase has been obtained as 49.3 kcal/mol at the UB3LYP/aug-cc-pVDZ level, while it has been computed as 45.7 kcal/mol at the UBLYP/aug-cc-pVDZ level (38,40,42) adopting the same molecular geometries. The energy distribution functions _QM(), and _QM() for the QM object in the solution and the reference systems are drawn in Fig. 4(a). It is shown that the energy distribution in the solution system locates on the lower region of the energy coordinate as compared with the reference system. It suggests that the electron affinity of the isoalloxazine ring is enhanced by the hydration. The free energy contribution from the isoalloxazine ring has been estimated to be 39.5 kcal/mol. In Fig. 4(b) the energy distribution functions _MM(), and _MM() for the interactions between the excess charge and the water solvent are presented, in which the distributions _QM(), and _QM() are also shown. As expected the interaction between the excess charge and the water solvent is enhanced in the solution system as compared with the reference system. This is simply because the water molecules behave under the influence of the solute-solvent interaction in the solution system as differed from the reference system. What we have to stress here is that there exists a positive correlation between the two kinds of solute-solvent interactions. In other words, increase in the stabilization between the solute (excess charge) and water molecules results in the increase in the electron affinity of the isoalloxazine ring. This is a direct consequence of the fact that the electron drawing effect of the water molecules bonded at the two oxygen atoms of the isoalloxazine ring enhances the electron affinity of the ring. It should be noted that the same effect can be expected for the isoalloxazine ring bound to the apoprotein. As shown in Fig. 5 two amino residues, Met122 and Phe487, are found to form distinctive hydrogen bonds at the two oxygen atoms. The reduction free energy due to the hydration has been computed as 40.6 kcal/mol which includes the free energy contribution of 10.9 kcal/mol estimated by the Born's equation (50) from bulk water outside the spherical cavity. Thus, the total reduction free energy has been obtained as 80.1 kcal/mol by the sum of these contributions.

4.3. Reduction free energy of FAD embedded in a protein system

It is straightforward to extend the method which identifies the excess electron as a solute to the computation of reduction free energy _red of the cholesterol oxidase 1B4V. The major difference that arises when the system is changed from the water solution to the protein system is that the number of kinds of solvents increases from two to four. The solvents in the protein system are, specifically, apoprotein that encloses the FAD, water molecules surrounding the protein, the isoalloxazine ring to be describe by Kohn-Sham DFT, and the part of FAD excluding the isoalloxazine ring. Correspondingly, we construct the energy distribution functions for the excess charge, _PRO(), _WAT(), _QM(), and _FAD(). The total energy distribution function () is simply given by the sum of these distributions. Owing to the fact that the energy distribution function can be divided into components, the reduction free energy can also be decomposed formally into contributions from various solvents. As noted in the previous subsection, the solution system is identified as the system in which the configurations of the solvent molecules are sampled under the influence of the excess charge located on the isoalloxazine ring. In the reference system, on the contrary, solute-solvent interaction is being switched off in the molecular simulations to construct the energy distribution functions.

Here we describe the outline of the preliminary QM/MM-ER simulations for the protein system. The whole of the isoalloxazine ring is described quantum chemically, and the remaining part of FAD as well as the apoprotein is represented by the CHARMM 22 force field. The boundary of the QM and the MM subsystem is described by a link atom. The apoprotein is immersed in a water droplet containing 8869 TIP3P molecules which are confined within a spherical cavity with a van der Waals wall. The center of the cavity is set at the center of mass of the isoalloxazine ring. In Fig. 6 the whole of the system is illustrated. The total charge of the protein is +3 in the unit of the elementary charge for the oxidized system and it is +2 for the reduced state.

The energy distribution functions for the solution and the reference systems are shown in Fig. 7. It can be readily recognized that there are four peaks in each distribution corresponding to the fact that there are four kinds of solvent. The free energy contributions from these solvents have been computed in terms of the individual energy distributions. It turns out that the largest contribution to the reduction free energy _red is given by the apoprotein (149.7 kcal/mol). The major origin for such the large stabilization of the excess charge can be attributed to the charge state of the apoprotein of which total charge amounts to +5 in the unit of the elementary charge. The secondary largest contribution is the free energy _red(WAT) due to the water molecules and it has been computed as 57.5 kcal/mol. It is surprising that the absolute value of free energy _red(WAT) in the protein system is larger than that in the water solution (40.6 kcal/mol), although it can be naturally speculated that substantial amount of water will be excluded from the vicinity of isoalloxazine ring by the protein. A possible fascinating interpretation for this phenomenon is that the protein induces the favorable hydration for the excess charge. A more detailed analysis must be done to clarify the underlying mechanism in the near future. The free energy changes due to the isoalloxazine ring and FAD without the ring have been computed as 31.9 kcal/mol and 55.1 kcal/mol, respectively. By utilizing the Born's equation we have approximately computed the free energy contribution _red(Born) from the bulk water. We have estimated _red(Born) as 21.0 kcal/mol for the charge variation of the system from +3 to +2. Thus, the free energy can be decomposed into terms corresponding to the individual energy distribution functions for the various kinds of solvent. We should note, however, that such the decomposition is possible only in the formulation since the effects of the many-body correlation among the solvents are inseparably involved in () in Eq. (3.23). Anyway, as a sum of these free energy components, the total reduction free energy has been obtained as _red = 163.1 kcal/mol by the QM/MM-ER approach, which is much larger than that computed in the aqueous solution (80.1 kcal/mol). The large difference in the free energies can be readily understood by attributing the origin to the large positive charge of the protein. According to the experimental measurement, however, the reduction free energy in the 1B4V system was estimated to be 90.6 kcal/mol (51). In our simulations the large positive charge of the system has not been compensated by the counter ions. On the other hand, it is likely that the charge of the protein in the experiment is being neutralized by counter ions coming from the substance added as a buffer for instance. Clearly it is desirable to compute the reduction free energy involving the counter ions as a variation to investigate the source of the deviation between the present computation and the experiment. However, the problem how to treat the charge of the protein in the simulation seems to be still unresolved to the best of author's knowledge. We have applied the QM/MM-ER method to compute the redox potential of cofactor embedded in an apoprotein. One of the notable features of the present approach is that the ensemble average is needed only for the initial and the final states of the reduction process in contrast to the free energy perturbation scheme owing to the use of the approximate functional based on the density functional theory in the energy representation. Moreover, it should be stressed that the convergence in the energy distribution function is very fast as compared with the spatial distribution functions, which makes a substantial advantage over the methods of the quantum chemical simulation in the free energy computation.

5. PERSPECTIVE

The free energy change plays an essential role in the exploration of the reaction mechanism in biological system since it governs the reaction pathway. Therefore, it is quite desirable to develop an efficient algorithm of quantum chemistry in combination with a statistical mechanical theory. Our recent development to solve the problem combines the hybrid quantum mechanical/molecular mechanical method with the theory of energy representation (QM/MM-ER). A characteristic feature in the description of the QM subsystem is that Kohn-Sham DFT with the real-space grid approach is used to achieve high performance in the parallel computations. Due mainly to the local nature of the operators in the effective Hamiltonian in Kohn-Sham DFT, real-space grids is able to afford a high efficiency in the parallel computation executed on a cheap work station cluster with distributed memory architectures. The QM/MM method has been applied to various chemical reactions in aqueous solutions and the efficiency and the adequacy has been established. Further, in this article we have tried to generalize the method to biological systems. We have started the development with the application to the reduction process of cofactor FAD enclosed in protein. For this purpose, we have taken a new approach where only the excess charge to be attached on the FAD is treated as a solute in the framework of the method of energy representation. Such a treatment is advantageous when it is applied to protein system because it enables ones to circumvent the direct evaluation of the protein's solvation free energy which may amount to a large negative value. The adequacy of the method has been examined for the reaction in aqueous solution by performing the conventional simulation based on the thermodynamic cycle. The reduction free energy of the active site of FAD, isoalloxazine ring, immersed in ambient water has been computed as 80.1 kcal/mol by the novel approach. Then, we have applied the method to cholesterol oxidase, 1B4V which contains FAD. The reduction free energy of FAD in 1B4V has been obtained as 163.1 kcal/mol. However, it was estimated to be 90.6 kcal/mol according to an experimental measurement. In this stage, we attribute such a serious deviation to the treatment of the large positive charge of the protein in the present simulation where the existence of the counter ions to compensate the positivity has not been assumed.

The approach presented in this review can be generally applied to the redox reaction in biological system, however, its application is limited only to the redox reaction. It is obvious that another procedure must be provided to compute free energy change for another sort of chemical reaction occurring in biological system. We categorize various kinds of chemical reactions into groups in order to adopt the method of QM/MM-ER. We are now going to construct a unified methodology for each group of chemical reactions. Our next subject to be studied is the proton affinity of the side chain in an amino residue embedded in a protein. pKa of the amino residue is obviously an important information which plays an essential role in determining the reaction mechanism and the protein structure as well. The present approach can also be applied to the computation of the ligand-binding free energy at the active site of a protein. The binding free energy is considered to be one of the most important quantity in the rational drag design. It is obvious that development of an efficient algorithm to determine the substrate-protein binding free energy gives a great impact in computational screening of the candidates. Thus, it is expected that the application of the QM/MM-ER approach will give substantial contributions to wide variety of fundamental chemical events in biological systems.

6. ACKNOWLEDGEMENT

This work was supported by Grant-in-Aid for Scientific Research (No. 18031022). It was also supported by the Next Generation Super Computing Project, Nanoscience Program, Mext, Japan. H. T. is very grateful to Prof. N. Matubayasi in Kyoto university for his valuable suggestions about the method of the energy representation and fruitful discussions. H. T. is also grateful to Prof. M. Nakano in Osaka university for his continuous encouragement to this work.

7. REFERENCES

1. A. Szabo and N. S. Ostlund: Modern Quantum Chemistry, Macmillan, New York. (1982)

2. R. G. Parr and W. Yang: Density Functional Theory of Atoms and Molecules, Oxford University Press, Oxford. (1989)

3. P. Hohenberg and W. Kohn: Inhomogeneous electron gas. Phys Rev 136, B864-B871 (1964)
doi:10.1103/PhysRev.136.B864

4. W. Kohn and L. J. Sham: Self-consistent equations including exchange and correlation effects. Phys Rev 140, A1133-A1138 (1965)
doi:10.1103/PhysRev.140.A1133

5. R. Car and M. Parrinello: Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys Rev Lett 55, 2471-2474 (1985)
doi:10.1103/PhysRevLett.55.2471

6. J. Gao and M. A. Thompson Ed.: Combined Quantum Mechanical and Molecular Mechanical Methods, American Chemical Society, Washington D.C. (1998)

7. M. F. Ruiz-Lopez Ed.: Combined QM/MM Calculations in Chemistry and Biochemistry, J Mol Struct (Theochem), special issue. 632, 1-336 (2003)

8. H. Takahashi, T. Hori, H. Hashimoto, and T. Nitta: A hybrid QM/MM method employing real space grids for QM water in the TIP4P water solvents. J Comput Chem 22, 1252-1261 (2001)
doi:10.1002/jcc.1082

9. T. Hori, H. Takahashi, and T. Nitta: Hybrid QM/MM molecular dynamics simulations for an ionic SN2 reaction in the supercritical water: OH^- + CH₃Cl CH₃OH + Cl^-. J Comput Chem 24, 209-221 (2003)
doi:10.1002/jcc.10134

10. H. Takahashi, S. Takei, T. Hori, and T. Nitta: A real-space-grid QM/MM study on the ionic/radical association reaction in aqueous phase: HCHO+OH→HCHO-OH. J Mol Struct (Theochem), 632, 185-195 (2003)
doi:10.1016/S0166-1280(03)00298-7

11. H. Takahashi, H. Hashimoto, and T. Nitta: Quantum mechanical/molecular mechanical studies of a novel reaction catalyzed by proton transfers in ambient and supercritical states of water. J Chem Phys 119, 7964-7971 (2003)
doi:10.1063/1.1610440

12. T. Hori, H. Takahashi, and T. Nitta: Hybrid quantum chemical studies for the methanol formation reaction assisted by the proton transfer mechanism in supercritical water: CH3Cl + nH2OCH3OH + HCl + (n-1)H2O. J Chem Phys 119, 8492-8499 (2003)
doi:10.1063/1.1611175

13. T. L. Beck: Real-space mesh techniques in density-functional theory. Rev Mod Phys 72, 1041-1080 (2000)
doi:10.1103/RevModPhys.72.1041

14. K. Hirose, T. Ono, Y. Fujimoto, and S. Tsukamoto: First-Principles Calculations in Real-Space Formalism, Imperial College Press, London. (2005)

15. J. P. Hansen and I. R. Macdonald: Theory of Simple Liquids, 2nd ed. Academic Press, London. (1986)

16. M. P. Allen and D. J. Tildesley: Computer Simulation of Liquids, Oxford University Press, Oxford. (1987)

17. D. Frenkel and B. Smit: Understanding Molecular Simulation, Academic Press, London. (2002)

18. J. Tomasi, B. Mennucci and R. Cammi: Quantum Mechanical Continuum Solvation Models. Chem Rev 105, 2999-3093 (2005)
doi:10.1021/cr9904009

19.D. Chandler and H. C. Andersen: Optimized Cluster Expansions for Classical Fluids. II. Theory of Molecular Liquids. J Chem Phys 57, 1930-1937 (1972)
doi:10.1063/1.1678513

20. F. Hirata and P. J. Rossky: An extended RISM equation for molecular polar fluids. Chem Phys Lett 83, 329-334 (1981)
doi:10.1016/0009-2614(81)85474-7

21. S. Ten-no, F. Hirata, and S. Kato: A hybrid approach for the solvent effect on the electronic structure of a solute based on the RISM and Hartree-Fock equations. Chem Phys Lett 214, 391-396 (1993)
doi:10.1016/0009-2614(93)85655-8

22. H. Takahashi, N. Matubayasi, M. Nakahara, and T. Nitta: A quantum chemical approach to the free energy calculations in condensed systems: The QM/MM method combined with the theory of energy representation. J Chem Phys 121, 3989-3999 (2004)
doi:10.1063/1.1774981

23. N. Matubayasi, and M. Nakahara: Theory of solutions in the energetic representation. I. Formulation. J Chem Phys 113, 6070-6081 (2000)
doi:10.1063/1.1309013

24.N. Matubayasi, and M. Nakahara: Theory of solutions in the energy representation. II. Functional for the chemical potential. J Chem Phys 117, 3605-3616 (2002)
doi:10.1063/1.1495850

N. Matubayasi, and M. Nakahara: Theory of solutions in the energy representation. II. Functional for the chemical potential. J Chem Phys 118, 2446 (2003)
doi:10.1063/1.1533752

25. N. Matubayasi, and M. Nakahara: Theory of solutions in the energy representation. III. Treatment of the molecular flexibility. J Chem Phys 119, 9686-9702 (2003)
doi:10.1063/1.1613938

26. H. Takahashi, W. Satou, T. Hori, and T. Nitta: An application of the novel quantum mechanical/molecular mechanical method combined with the theory of energy representation: An ionic dissociation of a water molecule in the supercritical water. J Chem Phys 122, 044504-044512 (2005)
doi:10.1063/1.1839858

27. H. Takahashi, Y. Kawashima, and T. Nitta: A novel quantum mechanical/molecular mechanical approach to the free energy calculation for isomerization of glycine in aqueous solution. J Chem Phys 123, 124504-124512 (2005)
doi:10.1063/1.2008234

28. T. Hori, H. Takahashi, M. Nakano, T. Nitta, and W. Yang: A QM/MM study combined with the theory of energy representation: Solvation free energies for anti/syn acetic acids in aqueous solution. Chem Phys Lett 419, 240-244 (2006)
doi:10.1016/j.cplett.2005.11.096

29. T. Hori, H. Takahashi, S. Furukawa, M. Nakano, and W. Yang: Computational Study on the Relative Acidity of Acetic Acid by the QM/MM Method Combined with the Theory of Energy Representation. J Phys Chem B 111,581-588 (2007)
doi:10.1021/jp066334d

30. H. Takahashi, K. Tanabe, M. Aketa, R. Kishi, S. Furukawa, and M. Nakano: Novel quantum mechanical/molecular mechanical method combined with the theory of energy representation: Free energy calculation for the Beckmann rearrangement promoted by proton transfers in the supercritical water. J Chem Phys 126, 084508 (2007)
doi:10.1063/1.2566834

31.J. R. Chelikowsky, N. Troullier, and Y. Saad: Finite-difference-pseudopotential method: Electronic structure calculations without a basis. Phys Rev Lett 72, 1240-1243 (1994)
doi:10.1103/PhysRevLett.72.1240

32. J. R. Chelikowsky, N. Troullier, K.Wu, and Y. Saad: Higher-order finite-difference pseudopotential method: An application to diatomic molecules. Phys. Rev. B 50, 11355-11364 (1994)
doi:10.1103/PhysRevB.50.11355

33. X. Jing, N. Troullier, D. Dean, N. Binggeli, J. R. Chelikowsky, K. Wu, and Y. Saad: Ab initio molecular-dynamics simulations of Si clusters using the higher-order finite-difference-pseudopotential method. Phys Rev B 50, 12234-12237 (1994)
doi:10.1103/PhysRevB.50.12234

34. A. D. Becke: A multicenter numerical integration scheme for polyatomic molecules. J Chem Phys 88, 2547-2553 (1988)
doi:10.1063/1.454033

35. G. B. Bachelet, D. R. Haman, and M. Schlüter: Pseudopotentials that work: From H to Pu. Phys Rev B 26, 4199-4228 (1982)
doi:10.1103/PhysRevB.26.4199

36. L. Kleinman and D. M. Bylander: Efficacious form for model pseudopotentials. Phys Rev Lett 48, 1425-1428 (1982)
doi:10.1103/PhysRevLett.48.1425

37. T. Ono and K. Hirose: Timesaving double-grid method for real-space electronic-structure calculations. Phys Rev Lett 82, 5016-5019 (1999)
doi:10.1103/PhysRevLett.82.5016

38. T. H. Dunning, Jr: Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J Chem Phys 90, 1007-1023 (1989)
doi:10.1063/1.456153

39. P. A. M. Dirac: Notes on exchange phenomena in the Thomas atom. Proc Cambridge Phil Soc 26, 376-385 (1930)

40. A. D. Becke: Density-functional exchange-energy approximation with correct asymptotic behavior. Phys Rev A 38, 3098-3100 (1988)
doi:10.1103/PhysRevA.38.3098

41. A. D. Becke: Density-functional thermochemistry. III. The role of exact exchange. J Chem Phys 98, 5648-5652 (1993)
doi:10.1063/1.464913

42. C. Lee, W. Yang, and R. G. Parr: Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys Rev B 37, 785-789 (1988)
doi:10.1103/PhysRevB.37.785

43. D. Voet, and J. G. Voet: Biochemistry 3rd ed., John Wiley & Sons, Inc., New Jersey. (2004)

44. H. Takahashi, H. Ohno, R. Kishi, M. Nakano, N. Matubayasi: Computation of the reduction free Energy of coenzyme in aqueous solution by the QM/MM-ER method. Chem Phys Lett to be published (2008)
doi:10.1016/j.cplett.2008.03.038

45. L. Guohui, X. Zhang, and Q. Cui: Free energy perturbation calculations with combined QM/MM potentials complications, simplifications, and applications to redox potential calculations. J Phys Chem B 107, 8643-8653 (2003)
doi:10.1021/jp034286g

46. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seigert: Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys Rev B 58 7260-7268 (1998)
doi:10.1103/PhysRevB.58.7260

47. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery. Jr., T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez, J. A. Pople: Gaussian 03, Revision B.05, Gaussian, Inc., Pittsburgh PA, 2003.

48. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein: Comparison of simple potential functions for simulating liquid water. J Chem Phys 79, 926-935 (1983)
doi:10.1063/1.445869

49. A. D. MacKerell, J. Wiorkiewicz-Kuczera, and M. Karplus: An all-atom empirical energy function for the simulation of nucleic acids. J Am Chem Soc 117, 11946-11975 (1995)
doi:10.1021/ja00153a017

A. D. MacKerell, Jr., D. Bashford, M. Bellott, R. L. Dunbrack, Jr., J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, III, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, and M. Karplus: All-atom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B, 102, 3586-3616 (1998)
doi:10.1021/jp973084f

50. M. Born: Volumen und hydratationswärme der ionen. Z Phys 1, 45-48 (1920)
doi:10.1007/BF01881023

51. Y. Yin, N. S. Sampson, A. Vrielink, P. I. Lario: The presence of a hydrogen bond between Asparagine 485 and the π system of FAD modulates the redox potential in the reaction catalyzed by cholesterol oxidase. Biochem 40, 13779-13787 (2001)
doi:10.1021/bi010843i

Abbreviations: ATP: adenosine triphosphate, aug-cc-pVDZ: correlation-consistent-valence-double-zeta basis set with polarization augmented by diffuse functions, BLYP: Becke, Lee, Yang, and Parr, B3LYP: Becke's three parameter model with LYP functional, CG: conjugate gradient, CI: configuration interaction, CP: Car-Parrinello, DFT: density functional theory, ER: energy representation, FAD: Flavin Adenine Dinucleotide, FEP: free energy perturbation, FFT: fast Fourier transformation, GGA: generalized gradient approximation, HF: Hartree-Fock, HNC: hypernetted chain, LCAO: linear combination of atomic orbitals, LDA: local density approximations, LJ: Lennard- Jones, Met: methionine, MO: molecular orbitals, PCM: polarizable continuum model, Phe: phenylalanine, PY: Percus-Yevick, QM/MM: quantum mechanical / molecular mechanical method, RISM: reference interaction site model, SCC-DFTB: self-consistent-charge density-functional-tight-binding, SCF: self-consistent field, SD: steepest descent, TI: thermodynamic integration, WAT: water

Key Words: Free energy, Quantum chemistry, Kohn-Sham DFT, QM/MM, Theory of solutions, Theory of energy representation , Review