Efficient Computational Research Protocol to ... - ACS Publications

Article pubs.acs.org/JPCB

Efficient Computational Research Protocol to Survey Free Energy Surface for Solution Chemical Reaction in the QM/MM Framework: The FEG-ER Methodology and Its Application to Isomerization Reaction of Glycine in Aqueous Solution Norio Takenaka,†,‡ Yukichi Kitamura,†,§ and Masataka Nagaoka*,†,‡,§ †

Graduate School of Information Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan ESICB, Kyoto University, Kyodai Katsura, Nishikyo-ku, Kyoto 615-8520, Japan § Core Research for Evolutional Science and Technology, Japan Science and Technology Agency, Honmachi, Kawaguchi 332-0012, Japan ‡

S Supporting Information *

ABSTRACT: In solution chemical reaction, we often need to consider a multidimensional free energy (FE) surface (FES) which is analogous to a Born− Oppenheimer potential energy surface. To survey the FES, an efficient computational research protocol is proposed within the QM/MM framework; (i) we first obtain some stable states (or transition states) involved by optimizing their structures on the FES, in a stepwise fashion, finally using the free energy gradient (FEG) method, and then (ii) we directly obtain the FE differences among any arbitrary states on the FES, efficiently by employing the QM/MM method with energy representation (ER), i.e., the QM/MM-ER method. To validate the calculation accuracy and efficiency, we applied the above FEG-ER methodology to a typical isomerization reaction of glycine in aqueous solution, and reproduced quite satisfactorily the experimental value of the reaction FE. Further, it was found that the structural relaxation of the solute in the QM/MM force field is not negligible to estimate correctly the FES. We believe that the present research protocol should become prevailing as one computational strategy and will play promising and important roles in solution chemistry toward solution reaction ergodography.

1. INTRODUCTION Solution chemical reaction always involves a variety of interactions among not only solute species but also a large number of solvent molecules. Quantum mechanical (QM) approaches can provide us the microscopically quite correct information regarding critical molecular structures, such as a stable state (SS) and transition state (TS). However, it still remains difficult to directly apply the QM approaches to large reaction systems such as solution systems and biomolecules. Thus, in the theoretical studies on such reaction systems, a quantum mechanical/molecular mechanical (QM/MM) method combined with molecular dynamics (MD) or Monte Carlo calculation is very useful to treat the whole systems in solution or enzymatic environment.1−3 In the gas phase, the chemical reactivity and reaction dynamics can be understood by certain characteristic features on a Born−Oppenheimer potential energy surface (PES) of such chemical reacting systems. However, because the PESs for the solution systems are too high-dimensional, it is almost impossible to obtain even SS of the solute molecule with respect to all degrees of freedom for the whole system. Under the circumstances, the free energy gradient (FEG) method with © 2016 American Chemical Society

the QM/MM-MD calculation (the QM/MM-FEG method) has been proposed to realize the full-atomic structural optimization of the solute molecule involved in solution.4−13 In the QM/MM-FEG method, the FEGs are utilized on the multidimensional free energy surface (FES) with respect to the structural coordinates of the QM solute for the free energy optimization by considering the contributions of a vast number of MM solvent molecules as the ensemble-averaged contribution under a certain thermodynamic environment. Moreover, because the Hessian on the FES, i.e., FE-Hessian, is similarly obtained in common with the FEG, a vibrational frequency analysis (VFA) with the QM/MM-MD calculation (QM/MMVFA) can be executed to obtain both effective normal modes and their vibrational frequencies in solution.9−13 Nevertheless, it is further difficult to obtain the free energy difference between two arbitrary states on the FES because its direct calculation is computationally demanding if such Special Issue: Bruce C. Garrett Festschrift Received: October 14, 2015 Revised: January 20, 2016 Published: January 21, 2016 2001

DOI: 10.1021/acs.jpcb.5b10061 J. Phys. Chem. B 2016, 120, 2001−2011

Article

The Journal of Physical Chemistry B

When investigating the electronic properties of an isolated molecule by QM calculations on the PES, one optimizes first the molecular structure and executes then some additional calculations necessary to estimate the specific physical chemical quantities by using the electron density analysis and VFA and so on. Moreover, a number of successive high-level theoretical methods are usually applied step by step if these methods should be required for such complicated systems that have difficulty in describing their subtle electronic states. It is true that such a general procedure should also be applicable to most solution reaction systems. In fact, as far as the same solute molecule is concerned, theoretical prediction of its properties of a solute molecule becomes usually more accurate in the following order: (i) the gas treatment, (ii) the DCM method, and (iii) the QM/MM method (high precision) if the same QM description for the solute is used in the three cases. Figure 1 shows the schematic flow diagram of the computational research protocol to study such physical and

standard techniques would be employed as free energy perturbation (FEP)14 and thermodynamic integration (TI)15 methods, which require estimation of free energies at many intermediate states between them. However, as for the free energy calculation at some fixed molecular structure, Matubayasi and Nakahara have proposed an energy representation (ER) method, reducing drastically the calculation cost, to compute the solvation free energy of the solute molecule, utilizing the MD calculations with an approximate functional for the free energy.16−18 Further, Takahashi and Matubayasi have developed the ER method combined with the QM/MMMD calculation (the QM/MM-ER method).19−22 On the basis of this method, one can calculate straightforwardly the solvation free energy of the QM solute as with the dielectric continuum model (DCM)23 without introducing the intermediate states. A recently proposed FEG method under a rigid body dynamics simulation, i.e., the FEG-RBD method,24 would be successfully applied with a benefit of the fast convergence of such costly FEGs. However, in the present article, we would like to propose another shortcut strategic protocol that, only if one could obtain some approximate reactant state (RS) and TS, or product state (PS), the activation or reaction free energy would be quite efficiently provided by employing the combined QM/MM-FEG and QM/MM-ER method, i.e., the QM/MMFEG-ER method, with the help of the QM/MM-VFA method. By utilizing such an efficient computational research protocol to survey FES for the solution chemical reaction, the “solution reaction ergodography”4 on FES would become, therefore, possible within the QM/MM framework beyond “reaction ergodography”25,26 on PES. We applied a QM/MM-FEG-ER method to an isomerization reaction of glycine in aqueous solution, which is a typical system to validate the calculation accuracy and efficiency of new computational methods.6,10,20,21,27−31 According to a previous quantum chemical study with the discrete continuum model,31 it was found that the addition of water molecules in the first solvation shell around the glycine molecule is required to reproduce the experimental reaction free energies of the tautomerization, showing that it is essential to include the explicit hydrogen-bonds (HBs) between glycine and water molecules. By using the present method, the reaction free energies can be accurately obtained by taking not only the HBs but also the thermodynamic fluctuation of water molecules under the finite temperature into consideration. In the present article, first, we investigated the geometries of the glycine molecule in solution through the structural optimization procedure. Then, we computed the solvation free energies of glycine, and estimated the reaction free energy for the isomerization reaction, discussing its dependence on the structural differences and theoretical calculation levels. After introducing generally the strategic methodology in section 2, we present calculation results with a few essential discussions in section 3, and draw conclusions including the future perspective in section 4.

Figure 1. Schematic flowchart of a computational scheme for the estimation of free energy in solution at any SS or TS on the FES.

chemical properties as solute equilibrium structures, the solute free energies, and vibrational frequencies at any SS or TS on the FES. At first, the structures of solute molecules are free energetically optimized in solution (FEG geometry optimization). In this procedure, one can obtain efficiently the stationary structures qSQM/MM on the FES and their potential energies EQM/MM by adopting those “final” structures optimized by the lower theoretical method as the initial guess geometry for the higher method. In particular, this is a natural but important step for the QM/MM-FEG method4−13 because it is required to execute a number of long QM/MM-MD calculations to collect many instantaneous forces at a certain molecular configurations in computing its FEG on the FES. After optimization of the solute structure via the QM/MMFEG method, the solvation free energy μ of the QM solute can be efficiently provided by employing the QM/MM-ER method.19−22 At the same time, because the structure of the QM solute is fixed in both the QM/MM-FEG and QM/MMER methods, the free energetic contribution from the intramolecular degrees of freedom for the solute as well as due to the quantum effect, i.e., zero point correction, is approximately obtained by using the effective vibrational frequencies {ν} in the QM/MM-VFA9−13 treatment, similarly

2. METHODS 2.1. QM/MM-FEG-ER Method. As an analogy to the energy calculations on the PES, we have proposed the QM/ MM-FEG-ER method that combines the QM/MM-FEG method4−13 with the QM/MM-ER method19−22 and is able to obtain the free energy changes between two arbitrary states on the FES, striking a good balance between the accuracy and calculation cost. 2002


Article


where the Helmholtz free energy A(qS) is a state function of not only thermodynamic variables (V, T) but also the solute structure qS. In eq 3, the brackets denote the equilibrium NVT ensemble average

as in the traditional treatment of the QM calculation for the gas molecule and in DCM.23 In the following sections (sections 2.2−2.4), the individual theoretical details of such three methods are shown by using the uniform variables among them. In the QM/MM-FEG-ER method, the free energy difference ΔAij can be obtained from one state i to another j on the FES as ΔAij = ΔEij + Δμij + ΔAijvib

⟨···⟩F =

(1)

(2a)

est vdW =⟨Ψ|Ĥ QM + Ĥ QM/MM|Ψ⟩ + Ĥ QM/MM + VMM

(2b)

=VSB + VMM

(2c)

ΔAi = Ai + 1 − Ai

S

F (q ) = −

∂A(q S) ∂q S

=−

(5b)

i

where qSi and qSi+1 are solute structures at the optimization steps i and i + 1, respectively. With application of the limitedmemory BFGS algorithm,33 the optimization cycle is repeated for the FEP treatment until the following condition, i.e., the zero gradient condition, is satisfied8 ∂VSB(q Si ) ∂q Si

≈0 q Si

(6)

2.3. QM/MM-VFA: Evaluation of Effective Vibration Frequencies on FES. The QM/MM-VFA has been realized by diagonalizing the mass-weighted “free energy (FE)” Hessian matrix in solution.9−13 For the analysis, the QM/MM-MD simulation is executed, calculating directly the elements of the FE Hessian matrix, which are described by the following equation. ∂ 2A(q S) = ∂q S∂q S

−

∂ 2VSB(q S) ∂q S∂q S

∂VSB(q S) ∂q S

qS

⎡ 1 ⎢ − kBT ⎢⎢ ⎣

∂VSB(q S) qS

∂q S

∂VSB(q S) ∂VSB(q S) ∂q S

⎤ ⎥ ⎥ qS ⎥ ⎦

∂q S

T

qS

T

(7)

Here the superscript T denotes the transposition. To reduce the simulation time, the FE Hessian matrix elements were calculated by using the following approximation ∂ 2A(q S) ∂q S∂q S

≈

∂ 2VSB(q S) ∂q S∂q S

qS

(8)

because the contribution of the second term in eq 7 was small enough in comparison to that of the first term.8−10 Then, one diagonalizes the “mass-weighted (mw)” FE Hessian matrix HFE ≡

∂ 2VSB(q S) ∂ mS q S∂ mS q S

≡ mS − 1/2 · qS

∂ 2A(q S) S − 1/2 ·m ∂q S∂q S (9a)

∂VSB(q S) ∂q S

(5a)

=−kBT ln⟨exp[−β {VSB(q Si + 1) − VSB(q Si )}]⟩q S

where Ĥ QM and Ĥ QM/MM stand for the standard Hamiltonian of the QM system and the interaction between the QM and the MM systems, respectively, and VMM is the potential energy of the MM system. The electrostatic interaction of the polarized QM charge density with an external charge distribution such as a set of MM point charges is estimated by the electronic embedding scheme32 which is implemented by adding the contribution of the MM point charges to the QM Hamiltonian, while the van der Waals interaction is typically described by a Lennard-Jones type potential. VSB denotes, therefore, the sum of both the potential energy of the solute and the interaction energy between the solute and solvent molecules. The QM/MM-FEG method4−13 has been developed for applications to identify the SS and TS structures and explore the FESs in solution reaction systems. In an equilibrium system, for example, the canonical (NVT) ensemble average of an observable is to be replaced as the time average calculated over a long enough NVT equilibrium QM/MM-MD simulation, provided the trajectory is ergodic. Then, the force on the NVT FES, FFE(qS), can be defined as a vector function of the positional coordinates of a solute structure qS and is equal to the time average of the instantaneous force vector acting on each atom of the NS atomic solute molecule as follows FE

(4)

where qB denotes the solvent positional coordinates as a whole, β is 1/kBT, and F is the arbitrary variable fixed during the simulation. By using the FEP method,14 the free energy difference at an optimization cycle i, ΔAi, is described as follows

where ΔEij is the potential energy change, Δμij is the change in the solvation free energy between the two states (see eq 13), and ΔAvib ij is the change due to the correction within the VFA, Avib (see eqs 10−12). If these two states correspond to an RS and a TS, or a PS in solution, respectively, such a difference would become the activation free energy, or the reaction free energy, respectively. The QM/MM-FEG method combined with the nudged elastic band (NEB) method, i.e., the QM/ MM-FEG-NEB method,9 therefore, should be a useful theoretical method in determining TSs for many chemically reacting systems on the FES, which are generally difficult to locate even in an isolated system. As an application example, we will show the capability of the method, focusing essentially on the reaction free energy of glycine isomerization in water. 2.2. QM/MM-FEG Method: Structural Optimization on FES. In the QM/MM treatment of reactive solute molecules in solution, it is common that only the reactive parts in the whole solution system are usually treated quantum mechanically, while the other part is molecular mechanically.1−3 The total system potential energy VT is, then, obtained by VT = ⟨Ψ|Ĥ QM + Ĥ QM/MM|Ψ⟩ + VMM

∫ dq B(···)exp(−βVT) ∫ dq B exp(−βVT)

qS

where mS denotes the following atomic mass matrix of the solute molecule

(3) 2003


Article

The Journal of Physical Chemistry B ⎛m S 0 ⎜ 1 ⎜ 0 mS 2 mS ≡ ⎜ ⎜⋮ ⋮ ⎜⎜ ⎝0 0

0 ⎞ ⎟ ⋯ 0 ⎟ ⎟ ⋱ 0 ⎟ ⎟ 0 m3SN S ⎟⎠

⎛ ρ e (ε ) ⎞ w e(ε) = −kBT ln⎜⎜ ue ⎟⎟ − ue(ε) ⎝ ρ0 (ε) ⎠

⋯

where u (ε) is the solute−solvent interaction (the direct part), ρeu(ε) and ρe0(ε) are the energy distribution functions in solution and the pure solvent system, respectively, and we(ε; λ) expresses the indirect part of the solute−solvent interaction (λ: a (set of) coupling parameter(s)). In eq 15a, the integration of we(ε; λ) is carried out with the combination of the Percus− Yevick (PY) and hypernetted chain (HNC) approximations.16−18 To estimate μf(qS) (in eq 13), we employed a recently proposed approach21 by Takahashi and Matubayasi, where the energy change H of the whole system due to the electron density polarization of the QM solute in the QM/MM method, is defined as

(9b)

and then one obtains the effective normal vibrational frequencies of the solute molecule in solution. Similarly, in the normal vibrational frequency analysis for a gas molecule, the intramolecular free energetic contributions Avib(qS) can be obtained by using a set of effective normal vibrational frequencies {νk(qS)|k = 1, ..., 3NS} at the solute structure qs, being defined under the second order approximation, as follows Avib(q S) = AZPE(q S) + AHA (q S)

(10)

H = ⟨Ψ|Ĥ QM|Ψ⟩ − ⟨Ψ0|Ĥ QM|Ψ0⟩

where AZPE is the contribution of the zero-point energy (ZPE) correction AZPE(q S) =

∑ k

1 hνk(q S) 2

+

k

Q (η) = ⟨δ(η − H )⟩q S

(12)

Q 0(η) = ⟨δ(η − H )⟩q S, n̅

μf =

(13)

dλw e(ε ; λ))(ρue (ε) − ρ0e (ε))]

⎞

⎞

⎝

0

⎠

(18a) (18b)

Here, we employed the weight function W(η) W (η) = C

21

Q (η)Q 0(η) Q (η) + Q 0(η)

(19)

where C is a normalization constant. In the present study, the QM/MM-ER method was implemented into the AMBER-Gaussian interface (AG-IF)34 by reference to the previous manuscripts on the ER method16−18 and QM/MM-ER method.19−22 In the implementation, the QM/MM-MD calculations with the fixed electron density n̅ were executed by loading the electron density matrix at any states in each MD step. Before such QM/ MM-MD calculations were performed, the electron density matrices at any states are needed by executing the molecular orbital (MO) or density functional theory (DFT) calculation using the Gaussian program.35 In the present QM/MM-ER implementation, several standard electronic state calculations are available, for example, the Hartree−Fock (HF), the second order Møller−Plesset perturbation theory (MP2), and a lot of DFT methods such as BLYP and B3LYP and so on. 2.5. Computational Details. We applied the QM/MMFEG-ER method to the isomerization reaction from a neutral

(14)

1

⎛

∫

∫ dε[(ρue (ε) − ρ0e (ε)) + βwe(ε)ρue (ε)

∫0

⎛

∫ ⎜⎜kBT ln⎜⎜⎝ QQ ((ηη)) ⎟⎟⎠ + η⎟⎟W (η) dη

= R (η)W (η) d η

where xj denotes the full coordinate of the jth solvent molecule, and v(n,̅ qS, xj) is the interaction potential energy between the solute with electron density n̅ at the solute structure qS and jth solvent molecule. The superscript e is attached to emphasize that a function is represented over the energy coordinate ε. Then, the solvation free energy μ̅ is estimated by

− (β

(17b)

respectively. In eq 17b, the average is taken over the samplings with both the fixed structure qs and electron density n.̅ Then, the solvation free energy μf is numerically obtained by

ρue (ε) = ⟨∑ δ(v(n ̅ , q S, x j) − ε)⟩u

μ ̅ = −kBT

(17a)

and

where μ̅ is regarded as the solvation free energy of a solute molecule with the electron density n.̅ On the other hand, μf is the contribution due to the electron density fluctuation around n̅ in solution. Then, the “energy” distribution function ρeu, i.e., the average number density distribution of those solvent molecules with a certain interaction energy, can be expressed as a function of interaction energy ε (as an independent coordinate) in the presence of a certain solute−solvent interaction u j

(16)

where Ψ0 is the electronic state wave function of an isolate solute molecule. Then, both the energy distribution function in solution Q(η) and that in a reference system Q0(η) as functions of the energy coordinate η are described as

(11)

2.4. QM/MM-ER Method: Computation of Solvation Free Energy of QM Solute. The QM/MM-ER method19−22 has been proposed by Takahashi and Matubayasi to compute the solvation free energy of the QM solute by using the QM/ MM-MD simulation. In the method, by introducing an electron density n̅ fixed at an arbitrary electronic state, the solvation free energy μ at a solute structure qS is obtained by μ(q S) = μ ̅ (q S) + μf (q S)

∑ (v(n , q S, x j) − v(n̅ , q S, x j)) j

while AHA is that within the harmonic approximation (HA) AHA (q S) = kBT ∑ ln{1 − exp( −hνk(q S)/kBT )}

(15b)

e

(15a)

with 2004


Article


Figure 2. Optimized structures of glycine in neutral form (NF) in the gas phase (a) and in aqueous solution (b and c) and in zwitterionic form (ZW) in aqueous solution (d and e), obtained by DCM (i.e., SCI-PCM method) and QM/MM-FEG method. Bond lengths are in Å, and bond angles are in degrees. The sign ± means that the solute glycine molecule takes two enantiomeric SS structures in aqueous solution.

field.42 The long-range electrostatic interactions between the MM molecules were calculated by the particle mesh Ewald method, and the cutoff distance of the LJ potential was 9.0 Å. The QM-MM electrostatic interactions were truncated by applying the cutoff method with the distance 12.0 Å. To execute the QM/MM-MD simulation, we used the AG-IF34 by combining the MD simulation program package Amber 9.043 with the ab initio QM calculation program Gaussian03.35 In the QM/MM-FEG calculation, at each FEG optimization cycle, a sampling run for the FEG averaging was executed for 30 ps with a time step of 1.0 fs. Using the optimized structure by the QM/MM-FEG method, to estimate the FE Hessian matrix for the QM/MM-VFA, the QM/MM-MD simulation was performed for 30 ps. On the other hand, in the QM/MM-ER calculation, to compute the solvation free energy μ̅, we executed 100 ps QM/MM-MD simulations with a fixed electron density n̅ to obtain the energy distribution function ρeu(ε) in the glycine aqueous solution system. At the same time, 200 ps MM-MD simulation was performed to obtain the energy distribution ρe0(ε) in a pure water system. In the computation of μf, the QM/MM-MD simulations were executed for 100 ps to obtain the energy distribution functions Q(η) and Q0(η).

form (NF) to a zwitterion (ZW) form of glycine molecule in aqueous solution, which is a typical example to test a number of new computational methods.6,10,20,21,27−31 According to the previous theoretical investigations,29−31 there are various conformers for both NF and ZW glycine molecules. Among them, we employed the cis-structure,10,36−38 which is often considered to correspond to RS and PS for the isomerization reaction of glycine molecule in aqueous solution.6,10,20,21,27−31 The QM/MM-MD simulations were performed for a system consisting of a single glycine molecule and 5306 water molecules in a cubic simulation cell with a side of 54.5 Å, which was adjusted initially by executing the classical MD simulation with NPT ensemble at 1 atm. In the classical MD simulation, we used the effective atomic charges of the solute glycine molecule obtained by the Merz−Kollman scheme39,40 at the B3LYP/6-31G(d) level of theory. The temperature was controlled to 300 K by the Berendsen algorithm. The threedimensional periodic boundary condition was imposed. For the QM/MM glycine−water solution system, we regarded a glycine molecule as the QM solute, which is described by the B3LYP/6-31G(d) level of theory. On the other hand, the whole solvent water molecules were taken as the MM part, and each water molecule was represented by the rigid TIP3P water model.41 These computational settings were determined on the basis of the treatment of a recent QM/MMMD simulation for the same chemical reacting system because the simulation was found to be able to reproduce reasonably the experimental value of the reaction free energy.29 The LJ interactions between the solute QM glycine and MM water molecules were obtained by using the general AMBER force

3. RESULTS AND DISCUSSION 3.1. Structural Optimization of Glycine in Aqueous Solution: DCM vs QM/MM-FEG Method. First, we discuss how the geometries of the solute glycine molecule change in aqueous solution. By using the present computational research protocol, the geometries of glycine molecules were optimized 2005


Article

The Journal of Physical Chemistry B in a stepwise fashion. Initially, a SCI-PCM method44 was used for the geometry optimization within the DCM treatment (see Figure 1). Figure 2 shows the optimized structures of both NF and ZW glycine molecules in the gas phase and those in aqueous solution by the DCM (i.e., SCI-PCM method) and the QM/MM-FEG method at the B3LYP/6-31G(d) level of theory, respectively. It should be noted that there is no SS structure of ZW glycine molecule in the gas phase because it is a matter of course that the ZW form should become stable only in aqueous solution. For those NF glycine molecules (Figure 2a−c), there was actually no symmetry in all the structures, even the Cs symmetry. However, those dihedral angles χ(C1−N4−O9− C7) were ±8.06° (gas phase), ± 2.13° (DCM), and ±2.06° (QM/MM), respectively, indicating that the NF glycine structure changes in solution so as to be closer to the planar one due to the solvent effect. On the other hand, in the ZW glycine molecule (Figure 2d,e), those structures took the Cs symmetric forms showing χ(C1−N4−O9-C7) = 0°, i.e., the molecular plane of the glycine molecule, as with the previous studies.7,27 In contrast to the NF glycine molecule, through the FEG optimization, both r(O9−H10)s and θ(C1−N4−H10)s increased significantly so as to increase the intramolecular distance between the −NH3 and −COO functional groups. Considering that the DCM underestimates the HB interactions between the solute and solvent molecules, such microscopic solvation must be essential for the structural optimization of the ionic species such as the ZW glycine molecule. Table 1 shows are the potential energies ΔE and dipole moments in the isolated state of the glycine molecule. Here, ΔE

solute glycine molecule can be more reasonably obtained by using the QM/MM-FEG method. 3.2. Computed Solvation Free Energy of Glycine in Aqueous Solution via QM/MM-ER Method Combined with QM/MM-FEG Method. We computed the solvation free energies of NF and ZW glycine molecules in aqueous solution by using the QM/MM-FEG-ER method. Specifically, we employed, as a solute structure for the QM/MM-ER method, the structure optimized by the QM/MM-FEG method. In Figure 3, the energy distribution functions ρeu(ε) and ρe0(ε) (see eq 14, etc.) for the (a) NF and (b) ZW glycine molecules were depicted and were obtained with use of two different electron densities of NF and ZW, such as that in the gas phase ngas (blue curves) and that by the SCI-PCM method nDCM (red curves), as the reference electron density n.̅ We have also shown the energy distribution functions ρen(ε) (black curves) obtained by calculating the solute−solvent interactions in aqueous solution with the “instantaneous” electron density n(t) during the 100 ps QM/MM-MD simulation. In the calculation, a pair interaction between the solute and a specific solvent molecule was calculated by utilizing the instantaneous electron density matrix obtained by the QM/MM calculation in each MD step while executing the QM/MM-MD simulation. Such a difference of ρeu(ε) from ρen(ε) indicates the deviation from the actual solute−solvent interactions in the QM/MM-MD simulation. In a comparison of ρeu(ε) with ρen(ε), the energy curves in the negative region became narrower in ρeu(ε) (see Figure 3a,b). This should be partly because (i) the polarization of the solute glycine is underestimated and (ii) the fluctuation of electron density is neglected. In fact, we found that the dipole moments of ngas (6.14 D, NF; 10.84 D, ZW) and nDCM (7.35 D, NF; 12.63 D, ZW) are underestimated in comparison to the average dipole moments during the QM/MM-MD simulation (8.35 ± 0.30 D, NF; 14.57 ± 0.26 D, ZW). It is natural that the pair interactions become stronger by using nDCM because the dipole moment of solute increases in the DCM. Accordingly, it is clearly understood that, for the ZW glycine molecule (Figure 3b), the energy peak position in the negative ε region of ρeu(ε) in use of nDCM shifted reasonably to a more negative direction of ε so as to become closer to that in ρen(ε). In Figure 4, the energy distribution functions Q(η) and Q0(η) (eqs 17a and 17b) were presented by using two such different electron densities as ngas (Figure 4a,c) and nDCM (Figure 4b,d). In the overlap region between Q(η) and Q0(η), the distribution functions W(η) and R(η) were estimated to compute μf (eqs 18a and 18b). As a result, it was found that the distribution functions Q(η) (blue curves) and Q0(η) (red curves) are clearly separated. This is because the statistical ensemble obtained by the QM/MM-MD simulation with the fixed electron density n̅ deviates considerably from that by the “straightforward” QM/MMMD simulation. In particular, when we employed ngas as the reference electron density, the overlap region between Q(η) and Q0(η) became rather small (compare Figure 4a,c to 4b,d). With the consideration that such separation is further enhanced in the ZW glycine molecule (Figure 4c), if such induced polarization of the solute might be quite large, it must be unsuitable to use ngas for the QM/MM-ER method. On the contrary, the use of nDCM was found to make the overlap between Q(η) and Q0(η) considerably improved, showing that R(η) becomes wider and contributes more to μf. In fact, it is true that the use of average electron density nave in the QM/ MM-MD simulation is most desirable.21 Nevertheless, in most

Table 1. Potential Energies ΔE’s (in kcal/mol) and Dipole Moments (in D) of NF and ZW Glycine Molecules in an Isolated Systema species

structure

ΔE

dipole moment

NF

gas DCM QM/MM-FEG DCM QM/MM-FEG

0 0.45 1.27 20.54 25.35

5.58 5.97 6.14 10.12 10.84

ZW a

All the values of potential energy are tabulated relative to that of NF glycine molecule in the gas phase.

values were evaluated relative to the potential energy of the NF glycine molecule in the gas phase. It was found that the dipole moments correlate positively with the potential energies. It can be understood that the electron density of glycine molecules distorts in aqueous solution so as to increase the dipole moment although the progress in the charge separation makes the potential energy unstable. It is understood that the increases in the dipole moment must lead to the stabilization of the free energy due to the enthalpy gain because the electrostatic interactions become stronger. Actually, it was reported that the FES of NF glycine conformers in aqueous solution becomes rather flat in comparison to the PES in the gas phase due to such propensity.38 By using the QM/MMFEG method, the dipole moments of the solute glycine molecule further increased because the increase in the dipole moment becomes more favorable under the explicit or atomistic solvent environment even if the potential energy is unstable. As a result, it is concluded that the geometry of the 2006


Article


Figure 3. Energy distribution functions for glycine in (a) neutral form (NF) and (b) zwitterionic form (ZW), obtained in the use of two such different electron densities as that in ngas (blue lines) and that by the SCI-PCM method nDCM (red lines). ρeu(ε) and ρe0 are the energy distribution functions in aqueous solution and pure water system, respectively. The energy distribution functions ρen(ε) in NF and ZW (black lines) obtained by using the instantaneous electron density are also shown. All the energy distribution functions are normalized by the bulk number density.

cases, it is a matter of course that the use of nDCM should be more useful instead of nave because nDCM can be easily obtained at the lower calculation cost than nave. Accordingly, we computed the solvation free energies μ as the sum of the principal contribution μ̅ and μf (see eqs 13, 15a, 15b, 18a, and 18b), summarized in Table 2. Compared with the principal solvation free energies μ̅ in-between ngas and nDCM, μ̅ with nDCM significantly increased in comparison to that with ngas due to the larger dipole moment of the solute glycine molecule. Nevertheless, it is very interesting to notice that, by adding μf to μ̅, both the total solvation free energies μ were obtained as the almost the same values, meaning that “the calculation results do not depend on the selection of the reference electron density n̅ as long as there is any overlap between Q(η) and Q0(η)”. This propensity is similar to that of the previous study by using the QM/MM-ER method.21 In addition, it can be said, therefore, that the QM/MM-ER method is compatible with the QM/ MM-FEG method, as far as the solvation free energy is concerned. 3.3. Reaction Free Energy for Isomerization Reaction of Glycine in Aqueous Solution: Accuracy and Efficiency of the Present Computational Research Protocol. According to the present computational research protocol (see eq 1), we evaluated the reaction free energies for the isomerization reaction from NF to ZW of glycine in aqueous solution (Table 3). In this section, we employed ngas as the reference electron density because it was found that the solvation free energies do not depend on the selection of electron density in the QM/MM-ER method as shown in Table 2. To discuss the importance of the structural relaxation by the QM/MM-FEG method, we also calculated the reaction free energy by using the solute glycine structure optimized by the DCM method (i.e., the SCI-PCM method). When we employed the optimized structure by DCM, ΔAvib was estimated by the VFA method in DCM. As a result, the reaction free energy in DCM was estimated to be 4.0 kcal/mol, which is extremely different from the experimental value, −7.3 or −7.7 kcal/mol.45,46 Such overestimation is the common

feature in the typical DCMs (see Table S1 in Supporting Information). On the other hand, by using the QM/MM-ER method, the reaction free energy became −5.8 kcal/mol and was significantly improved so as to be closer to the experimental values. Such difference is caused by the underestimation of the solvation free energy in DCM. This indicates that the explicit HB interaction, which is neglected in DCM, should be crucial to estimate the solvation free energy in aqueous solution. Further, by relaxing the molecular structure using the QM/ MM-FEG method, the reaction free energy was improved to become −7.8 kcal/mol due to the free energy stabilization of 2.0 kcal/mol. It is clear, therefore, that such an effect of the structural relaxation is not negligible to estimate quantitatively the reaction free energy in solution. In particular, this effect should become larger as the size of solute molecule increases. For comparison, in Table 3, the reaction free energy by using the FEP method with the QM/MM-MD calculation (i.e., the QM/MM-FEP method) is also shown (the simulation detail is described in Supporting Information). By using the QM/MMFEP method, the calculated reaction free energy became −7.5 kcal/mol, which is in very good agreement with that by the QM/MM-ER method (−7.8 kcal/mol). This clearly shows that the calculation accuracy of QM/MM-ER method is comparable with that of the QM/MM-FEP method. It is worth noting that, using the QM/MM-ER method, the calculation cost can be significantly reduced because the FEP treatment inevitably requires calculation of the ensemble averages for a number of intermediate states connecting between two SS structures. In the calculation, we have to consider an additional 80 intermediate states between NF and ZW of the glycine molecule. By a rough estimation, the present calculation time should become at least 1 order of magnitude less than that by directly using the QM/MM-FEP method. Namely, the QM/ MM-FEG-ER method can obtain the free energy change between two states such as RS and PS on the FES with much less computational cost, keeping its calculation accuracy. Next, we discuss the dependence on the selection of theoretical calculation levels to compute the reaction free 2007


Article


Figure 4. Energy distribution functions Q(η) (blue curves) and Q0(η) (red curves) in the aqueous solution and reference systems, respectively, of neutral form (NF) (a and b) and zwitterionic form (ZW) (c and d) of glycine in aqueous solution. The distribution functions R(η) and W(η) (eqs 18a and 19) are also shown as green curves and black ones (in the right axis), respectively.

Table 2. Solvation Free Energies μ and Their Components μ̅ and μf (in kcal/mol) of Glycine in Aqueous Solution Obtained by the QM/MM-ER Method species

density

NF

gas DCM gas DCM

ZW

μ̅ −10.5 −16.7 −34.7 −53.9

Table 3. Reaction Free Energies ΔA and Their Components ΔE, Δμ, and ΔAvib (in kcal/mol) from NF to ZW of Glycine in Aqueous Solutiona

μf

μ

method

structure

ΔE

Δμ

ΔAvib

ΔA

−5.3 1.0 −14.8 4.3

−15.8 −15.7 −49.5 −49.6

DCM QM/MM-ER QM/MM-ER QM/MM-FEP expt

DCM DCM QM/MM-FEG QM/MM-FEG

20.1 20.1 24.1

−18.3 −28.1 −33.7

2.2 2.2 1.8

4.0 −5.8 −7.8 −7.5b −7.3c, −7.7d

a

In the QM/MM-ER method, the solvation free energies were estimated by using ngas as the reference electron density. bComputational detail is described in Supporting Information. cExperimental value from ref 45. dExperimental value from ref 46.

energies. Table 4 shows reaction free energies ΔA and their components ΔE, Δμ, and ΔAvib by using three different calculation levels such as HF, MP2, and B3LYP with the 631G(d) basis set. As a result, the reaction free energies were estimated to be −10.5 and −14.9 kcal/mol at the HF and MP2 levels, respectively, which are overestimated in comparison to that at the B3LYP level (−7.8 kcal/mol). At the HF level, ΔE and Δμ (32.3, −44.3 kcal/mol) were significantly different from those at the B3LYP level (24.1, −33.7 kcal/mol). This is brought about by the large structural difference of the solute glycine molecule between the HF and B3LYP levels due to the

neglect of the electronic correlation effect (see Figure S2 in Supporting Information). On the other hand, while ΔE at the MP2 level (23.5 kcal/mol) was very similar to that at the B3LYP level (24.1 kcal/mol), Δμ at the former level was found to be significantly overestimated (−40.2 kcal/mol, MP2; −33.7 kcal/mol, B3LYP). In fact, the average dipole moments of the 2008


Article


multidimensional tunneling (VTST/MT)47 should be inevitable and should reinforce and generalize the present research protocol to be able to be applicable to any chemical reactions in aqueous solution,48,49 even in low temperature. In addition, in the present simulations, we did not consider the induced polarization of the MM solvent water molecules around the QM solute. It is considered, therefore, that the solvation free energy should be underestimated due to the neglect of the induced polarization of the water solvent. On the other hand, in the standard QM/MM-MD method which commonly considers only the solute molecule quantum mechanically, the QM-MM interactions should be artificially overestimated at the short distance due to the neglect of Pauli repulsion, inversely leading to the overestimation of the solvation free energy.50,51 Considering that these two errors could be canceled, the computed solvation free energy might sometimes become reasonable even if we neglected such effects. To directly resolve such a problem, it is straightforward that the closer MM solvent molecules around the QM solute should be included in widening the QM region because this problem is originating in the quantum-mechanical behaviors of a molecule. On the basis of such a strategy, we have been recently developing the number-adaptive multiscale (NAM) QM/MMMD method50,51 and the QM/MM-MD method combined with the fragment molecular orbital (FMO), i.e., the FMOQM/MM-MD method.52 However, even in the present computational research protocol, it is the structural optimization using the QM/MMFEG method that is most time-consuming under the current computational power. To further reduce the calculation cost, it should be reasonable, therefore, to utilize a number of excellent approximate methods developed so far, e.g., the ASEP/MD,53 ASEC,54 and mean-field QM/MM55 methods, which can reduce drastically the amount of QM calculations by utilizing the mean field approximation. On the other hand, it was recently reported that the calculation cost of the QM/MM-ER method is reduced using the second order perturbation theory (PT2).22 Moreover, even if the high-level theoretical methods might be required, the dual-level method that improves the free energies at the lower-level method with the correction at the high-level theory56,57 should be helpful for practical applications of the present methodology. Hence, it is no exaggeration to say that the present research protocol to survey FES must become realized generally in studying chemical reactions in the condensed state and the FES and its derivative quantities could be automatically explored by combination with these useful approximate methods, i.e., the solution reaction ergodography.

Table 4. Dependence on the Theoretical Calculation Levels To Compute the Reaction Free Energies ΔA and Their Components ΔE, Δμ, and ΔAvib (in kcal/mol)a level

ΔE

Δμ

ΔAvib

ΔA

HF MP2 B3LYP

32.3 23.5 24.1

−44.3 −40.2 −33.7

1.5 1.8 1.8

−10.5 −14.9 −7.8

a

The solvation free energies were estimated by using ngas as the reference electron density.

ZW glycine molecule during the QM/MM-MD simulations at the HF and MP2 levels were 15.76 ± 0.21 D and 14.93 ± 0.25 D, which are considerably overestimated in comparison to that at the B3LYP level (14.57 ± 0.25 D) if the 6-31G(d) basis set might be used in the three cases. It is considered, therefore, that the overestimation of Δμ at the HF and MP2 levels is mainly originating in the stronger electrostatic interactions between the solute glycine and solvent water molecules due to the larger dipole moments. However, the dependence on the calculation levels presented here might change if we would employ some larger basis sets. Nevertheless, it is true that the selection of the calculation level should strongly affect the reaction free energies. In particular, in the QM/MM method, because the numerical accuracy not only in the potential energy but also in the solvation free energy is each essential to estimate the reaction free energy, it is natural and necessary to determine more carefully the theoretical calculation level considering the balance between these two in solution.

4. CONCLUSIONS In the present study, we proposed an efficient computational research protocol, within the QM/MM method, to survey FES for the solution chemical reaction through the combination of the FEG, ER, and VFA methods. In the protocol, the FEG-ER methodology plays a central role to quantitatively estimate the reaction free energy and then was applied, as an example, to study the isomerization reaction of glycine in aqueous solution. As a result, the present research protocol using the QM/MMFEG-ER method was able to work accurately to estimate efficiently the reaction free energy without introducing the intermediate states like the FEP method. Further, it was clarified that the free energetic contribution of the structural relaxation in the QM/MM force field is not negligible to estimate correctly the reaction free energy. Nevertheless, it is a common matter of course that its calculation accuracy is strongly affected by the selection of theoretical calculation level, similar to the traditional electronic state calculations. Although we have shown here only one simple application in estimating the reaction free energy, the present research protocol can be straightforwardly applied also to evaluate the activation free energy. For example, only if any plausible TS structure on FES could be obtained by a theoretical method such as the QM/MM-FEG-NEB method,10 the activation free energy would be able to be estimated on the present computational strategy, as with the reaction free energy. However, if we would like to treat the proton transfer reaction such as the isomerization reaction of glycine in aqueous solution, the calculated reaction rate should be overestimated in comparison to the experimental one due to the neglect of nuclear quantum effects. To estimate these fine effects additionally, such a vibrational transition state theory including

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b10061. Figure S1 showing the free energy profile from NF to ZW glycine molecules with the QM/MM-FEP method and Figure S2 displaying optimized structures of glycine obtained by the HF and MP2 theoretical levels with the QM/MM-FEG method, Table S1 showing reaction free energies, and additional computational details (PDF) 2009


Article


■

Practical Aspects of Computational Chemisry III; Leszczynski, J., Shukla, M. K., Eds.; Springer: Heidelberg, 2014; pp 231−247. (13) Kitamura, Y.; Takenaka, N.; Koyano, Y.; Nagaoka, M. Free Energy Gradient Method and its Recent Related Developments: Free Energy Optimization and Vibrational Frequency Analysis in Solution. In Quantum Modeling of Complex Molecular Systems (Challenges and Advances in Computational Chemistry and Physics); Rivail, J.-L., RuizLopez, M., Assfeld, X., Eds.; Springer: Heidelberg, 2015; pp 219−252. (14) Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, 22, 1420−1426. (15) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3, 300−313. (16) Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energetic Representation. I. Formulation. J. Chem. Phys. 2000, 113, 6070−6081. (17) Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energy Representation. II. Functional for the Chemical Potential. J. Chem. Phys. 2002, 117, 3605−3616. (18) Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energy Representation. III. Treatment of the Molecular Flexibility. J. Chem. Phys. 2003, 119, 9686−9702. (19) Takahashi, H.; Matubayasi, N.; Nakahara, M.; Nitta, T. 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. 2004, 121, 3989− 3999. (20) Takahashi, H.; Kawashima, Y.; Nitta, T.; Matubayasi, N. A Novel Quantum Mechanical/Molecular Mechanical Approach to the Free Energy Calculation for Isomerization of Glycine in Aqueous Solution. J. Chem. Phys. 2005, 123, 124504. (21) Takahashi, H.; Omi, A.; Morita, A.; Matubayasi, N. Simple and Exact Approach to the Electronic Polarization Effect on the Solvation Free Energy: Formulation for Quantum-Mechanical/MolecularMechanical System and Its Applications to Aqueous Solutions. J. Chem. Phys. 2012, 136, 214503. (22) Suzuoka, D.; Takahashi, H.; Morita, A. Computation of the Free Energy Due to Electron Density Fluctuation of a Solute in Solution: A QM/MM Method with Perturbation Approach Combined with a Theory of Solutions. J. Chem. Phys. 2014, 140, 134111. (23) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3093. (24) Tao, P.; Sodt, A. J.; Shao, Y.; König, G.; Brooks, B. R. Computing the Free Energy along a Reaction Coordinate Using Rigid Body Dynamics. J. Chem. Theory Comput. 2014, 10, 4198−4207. (25) Fukui, K. The Path of Chemical Reactions - The IRC Approach. Acc. Chem. Res. 1981, 14, 363−368. (26) Fukui, K. The Role of Frontier Orbitals in Chemical Reactions. Angew. Chem., Int. Ed. Engl. 1982, 21, 801−809. (27) Tuñoń , I.; Silla, E.; Ruiz-López, M. F. On the Tautomerization Process of Glycine in Aqueous Solution. Chem. Phys. Lett. 2000, 321, 433−437. (28) Leung, K.; Rempe, S. B. Ab Initio Molecular Dynamics Study of Glycine Intramolecular Proton Transfer in Water. J. Chem. Phys. 2005, 122, 184506. (29) Lu, Z.; Zhamg, Y. Interfacing ab Initio Quantum Mechanical Method with Classical Drude Osillator Polarizable Model for Molecular Dynamics Simulation of Chemical Reactions. J. Chem. Theory Comput. 2008, 4, 1237−1248. (30) Yoshikawa, T.; Motegi, H.; Kakizaki, A.; Takayanagi, T.; Shiga, M.; Tachikawa, M. Path-Integral Molecular Dynamics Simulations of Glycine·(H2O)n (n = 1−7) Clusters on Semi-Empirical PM6 Potential Energy Surfaces. Chem. Phys. 2009, 365, 60−68. (31) Benbrahim, N.; Rahmouni, A.; Ruiz-López, M. F. A Theoretical Study of Medium Effects on the Structure of the Glycine Analogue Aminomethylphosphonic Acid. Phys. Chem. Chem. Phys. 2008, 10, 5624−5632.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone/fax: +81-52-7895623 URL: http://www.ncube.human.nagoya-u.ac.jp/. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by a Grant-in-Aid for Science Research from the Ministry of Education, Culture, Sport, Science and Technology (MEXT) in Japan; also by the Core Research for Evolutional Science and Technology (CREST) “Establishment of Molecular Technology towards the Creation of New Functions” of the Japan Science and Technology Agency (JST); and partially by the MEXT program “Elements Strategy Initiative to Form Core Research Center” (since 2012), Japan.

■

REFERENCES

(1) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. (2) Singh, U. C.; Kollman, P. A. A Combined ab Initio Quantum Mechanical and Molecular Mechanical Method for Carrying Out Simulations on Complex Molecular Systems: Applications to the CH3Cl + Cl− Exchange Reaction and Gas Phase Protonation of Polyethers. J. Comput. Chem. 1986, 7, 718−730. (3) Field, M. J.; Bash, P. A.; Karplus, M. A. A Combined Quantum Mechanical and Molecular Mechanical Potential for Molecular Dynamics Simulations. J. Comput. Chem. 1990, 11, 700−733. (4) Okuyama-Yoshida, N.; Nagaoka, M.; Yamabe, T. Transition-State Optimization on Free Energy Surfaces: Toward Solution Chemical Reaction Ergodography. Int. J. Quantum Chem. 1998, 70, 95−103. (5) Okuyama-Yoshida, N.; Nagaoka, M.; Yamabe, T. Potential Energy Function for Intramolecular Proton Transfer Reaction of Glycine in Aqueous Solution. J. Phys. Chem. A 1998, 102, 285−292. (6) Nagaoka, M.; Okuyama-Yoshida, N.; Yamabe. Origin of the Transition State on the Free Energy Surface: Intramolecular Proton Transfer Reaction of Glycine in Aqueous Solution. T. J. Phys. Chem. A 1998, 102, 8202−8208. (7) Okuyama-Yoshida, N.; Kataoka, K.; Nagaoka, M.; Yamabe, T. Structure Optimization via Free Energy Grandient Method: Application to Glycine Zwitterion in Aqueous Solution. J. Chem. Phys. 2000, 113, 3519−3524. (8) Nagaoka, M.; Nagae, Y.; Koyano, Y.; Oishi, Y. Transition-State Characterization of Ammonia Ionization Process in Aqueous Solution via Free-Energy Gradient Method. J. Phys. Chem. A 2006, 110, 4555− 4563. (9) Koyano, Y.; Takenaka, N.; Nakagawa, Y.; Nagaoka, M. On the Importance of Lennard-Jones Parameter Calibration in QM/MM Framework: Reaction Path Tracing via Free Energy Gradient Method for Ammonia Ionization Process in Aqueous Solution. Bull. Chem. Soc. Jpn. 2010, 83, 486−494. (10) Takenaka, N.; Kitamura, Y.; Koyano, Y.; Asada, T.; Nagaoka, M. Reaction Path Optimization and Vibrational Frequency Analysis via ab Initio QM/MM Free Energy Gradient (FEG) Method: Application to Isomerization Process of Glycine in Aqueous Solution. Theor. Chem. Acc. 2011, 130, 215−226. (11) Kitamura, Y.; Takenaka, N.; Koyano, Y.; Nagaoka, M. Dual Approach to Vibrational Spectra in Solution: Microscopic Influence of Hydrogen Bonding to the State of Motion of Glycine in Water. J. Chem. Theory Comput. 2014, 10, 3369−3378. (12) Georg, H. C.; Fernandes, T. S.; Canuto, S.; Takenaka, N.; Kitamura, Y.; Nagaoka, M. A Combination of the Sequential QM/MM and the Free Energy Gradient Methodologies with Applications. In 2010


Article

The Journal of Physical Chemistry B (32) Bakowies, D.; Thiel, W. Hybrid Models for Combined Quantum Mechanical and Molecular Mechanical Approaches. J. Phys. Chem. 1996, 100, 10580−10594. (33) Liu, D. C.; Nocedal, J. On the Limited Memory BFGS Method for Large Scale Optimization. Math. Prog. 1989, 45, 503−528. (34) Okamoto, T.; Yamada, K.; Koyano, Y.; Asada, T.; Koga, N.; Nagaoka, M. A Minimal Implementation of the AMBER−GAUSSIAN Interface for ab Initio QM/MM-MD Simulation. J. Comput. Chem. 2011, 32, 932−942. (35) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr; Vreven, T.; Kudin, K. N.; Burant, J. C.; et al. Gaussian03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (36) Császár, A. G. Conformers of Gaseous Glycine. J. Am. Chem. Soc. 1992, 114, 9568−9575. (37) Jensen, J. H.; Gordon, M. S. On the Number of Water Molecules Necessary To Stabilize the Glycine Zwitterion. J. Am. Chem. Soc. 1995, 117, 8159−8170. (38) Kitamura, Y.; Takenaka, N.; Koyano, Y.; Nagaoka, M. On the Smoothing of Free Energy Landscape of Solute Molecules in Solution: A Demonstration of the Stability of Glycine Conformers via Ab Initio QM/MM Free Energy Calculation. Chem. Phys. Lett. 2011, 514, 261− 266. (39) Singh, U. C.; Kollman, P. A. An Approach to Computing Electrostatic Charges for Molecules. J. Comput. Chem. 1984, 5, 129− 145. (40) Besler, B. H.; Merz, K. M., Jr; Kollman, P. A. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11, 431−439. (41) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; et al. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (42) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (43) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Pearlman, D. A.; Crowley, M.; et al. AMBER 9; University of California: San Francisco, 2006. (44) Foresman, J. B.; Keith, T. A.; Wiberg, K. B.; Snoonian, J.; Frisch, M. J. Solvent Effects 5. The Influence of Cavity Shape, Truncation of Electrostatics, and Electron Correlation on ab Initio Reaction Field Calculations. J. Phys. Chem. 1996, 100, 16098−16104. (45) Haberfield, P. What is the Energy Difference between H2NCH2CO2H and + H3NCH2CO2-? J. Chem. Educ. 1980, 57, 346−347. (46) Wada, G.; Tamura, E.; Okina, M.; Nakamura, M. On the Ratio of Zwitterion Form to Uncharged Form of Glycine at Equilibrium in Various Aqueous Media. Bull. Chem. Soc. Jpn. 1982, 55, 3064−3067. (47) Truhlar, D. G.; Isaacson, A. D.; Skodje, R. T.; Garrett, B. C. Incorporation of Quantum Effects in Generalized-Transition-State Theory. J. Phys. Chem. 1982, 86, 2252−2261. (48) Truhlar, D. G.; Liu, Y.-P.; Schenter, G. K.; Garrett, B. C. Tunneling in the Presence of a Bath: A Generalized Transition State Theory Approach. J. Phys. Chem. 1994, 98, 8396−8405. (49) McRae, R. P.; Schenter, G. K.; Garrett, B. C.; Svetlicic, Z.; Truhlar, D. G. Vibrational Transition State Theoy Evaluation of the Rate Constant for Proton Transfer in a Polar Solvent. J. Chem. Phys. 2001, 115, 8460−8480. (50) Takenaka, N.; Kitamura, Y.; Koyano, Y.; Nagaoka, M. The Number-Adaptive Multiscale QM/MM Molecular Dynamics Simulation: Application to Liquid Water. Chem. Phys. Lett. 2012, 524, 56− 61. (51) Takenaka, N.; Kitamura, Y.; Koyano, Y.; Nagaoka, M. An Improvement in Quantum Mechanical Description of Solute-Solvent Interactions in Condensed Systems via the Number-Adaptive Multiscale Quantum Mechanical/Molecular Mechanical-Molecular Dynamics Method: Application to Zwitterionic Glycine in Aqueous Solution. J. Chem. Phys. 2012, 137, 024501.

(52) Okamoto, T.; Ishikawa, T.; Koyano, Y.; Yamamoto, N.; Kuwata, K.; Nagaoka, M. A Minimal Implementation of the AMBER-PAICS Interface for ab Initio FMO-QM/MM-MD Simulation. Bull. Chem. Soc. Jpn. 2013, 86, 210−222. (53) Galván, I. F.; Sánchez, M. L.; Martín, M. E.; Olivares del Valle, F. J.; Aguilar, M. A. ASEP/MD: A Program for the Calculation of Solvent Effects Combining QM/MM Methods and the Mean Field Approximation. Comput. Phys. Commun. 2003, 155, 244−259. (54) Georg, H. C.; Canuto, S. Electronic Properties of Water in Liquid Environment. A Sequential QM/MM Study Using the Free Energy Gradient Method. J. Phys. Chem. B 2012, 116, 11247−11254. (55) Nakano, H.; Yamamoto, T. Variational Calculation of Quantum Mechanical/Molecular Mechanical Free Energy with Electronic Polarization of Solvent. J. Chem. Phys. 2012, 136, 134107. (56) Retegan, M.; Martins-Costa, M.; Ruiz-López, M. F. Free Energy Calculations using Dual-Level Born-Oppenheimer Molecular Dynamics. J. Chem. Phys. 2010, 133, 064103. (57) Martins-Costa, M. T. C.; Ruiz-López, M. F. Amino Acid Capture by Aqueous Interfaces. Implications for Biological Uptake. J. Phys. Chem. B 2013, 117, 12469−12474.

2011


Efficient Computational Research Protocol to ... - ACS Publications

Recommend Documents