^{1}and Takeshi Yamamoto

^{1,a)}

### Abstract

Quantum mechanical/molecular mechanical (QM/MM) free energy calculation presents a significant challenge due to an excessive number of QM calculations. A useful approach for reducing the computational cost is that based on the mean field approximation to the QM subsystem. Here, we describe such a mean-field QM/MM theory for electronically polarizable systems by starting from the Hartree product ansatz for the total system and invoking a variational principle of free energy. The MM part is then recast to a classical polarizable model by introducing the charge response kernel. Numerical test shows that the potential of mean force (PMF) thus obtained agrees quantitatively with that obtained from a direct QM/MM calculation, indicating the utility of self-consistent mean-field approximation. Next, we apply the obtained method to prototypical reactions in several qualitatively different solvents and make a systematic comparison of polarization effects. The results show that in aqueous solution the PMF does not depend very much on the water models employed, while in nonaqueous solutions the PMF is significantly affected by explicit polarization. For example, the free energy barrier for a phosphoryl dissociationreaction in acetone and cyclohexane is found to increase by more than 10 kcal/mol when switching the solvent model from an empirical to explicitly polarizable one. The reason for this is discussed based on the parametrization of empirical nonpolarizable models.

This work was supported by the Grant-in-Aid for Scientific Research (21350010) and that on Innovative Areas (21118508) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan.

I. INTRODUCTION

II. THEORY

A. QM/MM free energy for electronically polarizable systems

B. Mean field approximation to QM/MM free energy

III. COMPUTATIONAL DETAILS

A. Calculation of the PMF profile

B. Protocol of CRK-MD simulation

IV. RESULTS AND DISCUSSION

A. Finkelstein reaction

B. Menshutkin reaction

C. Phosphoryl dissociationreaction

V. CONCLUSIONS

### Key Topics

- Solvents
- 89.0
- Free energy
- 80.0
- Mean field theory
- 38.0
- Polarization
- 22.0
- Wave functions
- 21.0

## Figures

Potential energy profile of the Finkelstein reaction (Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−}) in the gas phase. Zero of energy is set at the infinitely separated reactants. The basis sets are 6-311++G(3*df*,3*p*) for CCSD(T) and MP2 and 6-31+G(*d*,*p*) for B3LYP, BHHLYP, and HF methods. For comparison, the energy profile at the HF/6-31G(*d*) level is also shown.

Potential energy profile of the Finkelstein reaction (Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−}) in the gas phase. Zero of energy is set at the infinitely separated reactants. The basis sets are 6-311++G(3*df*,3*p*) for CCSD(T) and MP2 and 6-31+G(*d*,*p*) for B3LYP, BHHLYP, and HF methods. For comparison, the energy profile at the HF/6-31G(*d*) level is also shown.

Potential of mean force (PMF) of the Finkelstein reaction Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−} in aqueous solution. The result obtained with the mean-field QM/MM method is compared to those obtained from direct QM/MM calculations (see text for details). The QM calculation is performed at the HF/6-31G(*d*) level. The water is represented by the TIP3P model. The circles depict raw data points obtained by integrating the free energy gradient in Eq. (30). The statistical error of the PMF is comparable to the width of the plotted curve (see the supplementary material^{116} for more details).

Potential of mean force (PMF) of the Finkelstein reaction Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−} in aqueous solution. The result obtained with the mean-field QM/MM method is compared to those obtained from direct QM/MM calculations (see text for details). The QM calculation is performed at the HF/6-31G(*d*) level. The water is represented by the TIP3P model. The circles depict raw data points obtained by integrating the free energy gradient in Eq. (30). The statistical error of the PMF is comparable to the width of the plotted curve (see the supplementary material^{116} for more details).

PMF of the Finkelstein reaction Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−} in water and acetone solutions calculated with the mean-field QM/MM method. Solvents are described with the TIP3P and OPLS models in the nonpolarizable case (dashed lines) and with the CRK model in the polarizable case (solid line). The QM calculation is performed at the BHHLYP/6-31+G(*d*,*p*) level.

PMF of the Finkelstein reaction Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−} in water and acetone solutions calculated with the mean-field QM/MM method. Solvents are described with the TIP3P and OPLS models in the nonpolarizable case (dashed lines) and with the CRK model in the polarizable case (solid line). The QM calculation is performed at the BHHLYP/6-31+G(*d*,*p*) level.

Solvent electrostatic potentials (ESPs) acting on the solute atoms, 〈*V* _{α}〉, for the Finkelstein reaction in (a) water and (b) acetone solutions. The results obtained with the polarizable and nonpolarizable solvent models are shown by solid and dashed lines, respectively. Cl* and Cl indicate the attacking and leaving chloride atoms, respectively. The curve labeled with CH_{3} shows the mean value of ESP acting on the methyl group. Panel (c) displays the solvation free energy calculated with the linear response approximation (see the main text). The ESP values in panels (a) and (b) include the contribution of the Wigner potential.

Solvent electrostatic potentials (ESPs) acting on the solute atoms, 〈*V* _{α}〉, for the Finkelstein reaction in (a) water and (b) acetone solutions. The results obtained with the polarizable and nonpolarizable solvent models are shown by solid and dashed lines, respectively. Cl* and Cl indicate the attacking and leaving chloride atoms, respectively. The curve labeled with CH_{3} shows the mean value of ESP acting on the methyl group. Panel (c) displays the solvation free energy calculated with the linear response approximation (see the main text). The ESP values in panels (a) and (b) include the contribution of the Wigner potential.

PMF of the Finkelstein reaction in various solvents calculated with (a) the mean-field QM/MM method and (b) the COSMO continuum solvation model. In panel (a), all the solvents are described with the polarizable CRK model. The QM calculation is performed at the BHHLYP/6-31G+(*d*,*p*) level.

PMF of the Finkelstein reaction in various solvents calculated with (a) the mean-field QM/MM method and (b) the COSMO continuum solvation model. In panel (a), all the solvents are described with the polarizable CRK model. The QM calculation is performed at the BHHLYP/6-31G+(*d*,*p*) level.

Potential energy profile of the Menshutkin reaction NH_{3} + CH_{3}Cl → NH_{3}CH_{3} ^{+} + Cl^{−} in the gas phase. The basis sets are 6-311++G(3*df*,3*p*) for CCSD(T) and MP2 methods and 6-31+G(*d*,*p*) otherwise.

Potential energy profile of the Menshutkin reaction NH_{3} + CH_{3}Cl → NH_{3}CH_{3} ^{+} + Cl^{−} in the gas phase. The basis sets are 6-311++G(3*df*,3*p*) for CCSD(T) and MP2 methods and 6-31+G(*d*,*p*) otherwise.

PMF of the Menshutkin reaction NH_{3} + CH_{3}Cl → NH_{3}CH_{3} ^{+} + Cl^{−} in water, DMF, and cyclohexane solutions calculated with the mean-field QM/MM method. Solvents are described with the TIP3P and OPLS models in the nonpolarizable case (dashed lines) and with the CRK model in the polarizable case (solid line). The QM calculation is performed at the BHHLYP/6-31+G(*d*,*p*) level.

PMF of the Menshutkin reaction NH_{3} + CH_{3}Cl → NH_{3}CH_{3} ^{+} + Cl^{−} in water, DMF, and cyclohexane solutions calculated with the mean-field QM/MM method. Solvents are described with the TIP3P and OPLS models in the nonpolarizable case (dashed lines) and with the CRK model in the polarizable case (solid line). The QM calculation is performed at the BHHLYP/6-31+G(*d*,*p*) level.

Solvent electrostatic potentials (ESPs) acting on the solute atoms, 〈*V* _{α}〉, for the Menshutkin reaction in (a) water, (b) DMF, and (c) cyclohexane solutions. The results obtained with the polarizable and nonpolarizable solvent models are shown by solid and dashed lines, respectively. Panel (d) displays the solvation free energy calculated with the linear response approximation (see the main text). In panel (c), the ESP value calculated with the OPLS-UA model for cyclohexane is identically zero because the latter has no partial charge on the united CH_{2} atoms.

Solvent electrostatic potentials (ESPs) acting on the solute atoms, 〈*V* _{α}〉, for the Menshutkin reaction in (a) water, (b) DMF, and (c) cyclohexane solutions. The results obtained with the polarizable and nonpolarizable solvent models are shown by solid and dashed lines, respectively. Panel (d) displays the solvation free energy calculated with the linear response approximation (see the main text). In panel (c), the ESP value calculated with the OPLS-UA model for cyclohexane is identically zero because the latter has no partial charge on the united CH_{2} atoms.

PMF of the Menshutkin reaction in various solvents calculated with (a) the mean-field QM/MM method and (b) the COSMO continuum solvation model. In panel (a), all the solvents are described with the polarizable CRK model. The QM calculation is performed at the BHHLYP/6-31+G(*d*,*p*) level.

PMF of the Menshutkin reaction in various solvents calculated with (a) the mean-field QM/MM method and (b) the COSMO continuum solvation model. In panel (a), all the solvents are described with the polarizable CRK model. The QM calculation is performed at the BHHLYP/6-31+G(*d*,*p*) level.

Potential energy profiles of the phosphoryl dissociation reaction CH_{3}OPO_{3} ^{2 −} → CH_{3}O^{−} + PO_{3} ^{−} in the gas phase. The basis sets are 6-311++G(3*df*,3*p*) for CCSD(T) and MP2 methods and 6-31+G(*d*,*p*) otherwise.

Potential energy profiles of the phosphoryl dissociation reaction CH_{3}OPO_{3} ^{2 −} → CH_{3}O^{−} + PO_{3} ^{−} in the gas phase. The basis sets are 6-311++G(3*df*,3*p*) for CCSD(T) and MP2 methods and 6-31+G(*d*,*p*) otherwise.

PMF of the phosphoryl dissociation reaction CH_{3}OPO_{3} ^{2 −} → CH_{3}O^{−} + PO_{3} ^{−} in (a) acetone and (b) cyclohexane solutions calculated with the mean-field QM/MM method. For comparison, PMF obtained with the COSMO model and the potential energy profile in the gas phase are also shown. The QM calculation is performed at the MPW1PW91/6-31+G(*d*,*p*) level.

PMF of the phosphoryl dissociation reaction CH_{3}OPO_{3} ^{2 −} → CH_{3}O^{−} + PO_{3} ^{−} in (a) acetone and (b) cyclohexane solutions calculated with the mean-field QM/MM method. For comparison, PMF obtained with the COSMO model and the potential energy profile in the gas phase are also shown. The QM calculation is performed at the MPW1PW91/6-31+G(*d*,*p*) level.

Solvent electrostatic potentials (ESPs) acting on the solute atoms, 〈*V* _{α}〉, for the phosphoryl dissociation reaction in (a) acetone and (b) cyclohexane solutions. The results obtained with the polarizable and nonpolarizable solvent models are shown by solid and dashed lines, respectively. Panel (c) displays the solvation free energy calculated with the linear response approximation (see the main text). The ESP values in panels (a) and (b) include the contribution of the Wigner potential.

Solvent electrostatic potentials (ESPs) acting on the solute atoms, 〈*V* _{α}〉, for the phosphoryl dissociation reaction in (a) acetone and (b) cyclohexane solutions. The results obtained with the polarizable and nonpolarizable solvent models are shown by solid and dashed lines, respectively. Panel (c) displays the solvation free energy calculated with the linear response approximation (see the main text). The ESP values in panels (a) and (b) include the contribution of the Wigner potential.

## Tables

Free energy barrier Δ*A* ^{‡} (in kcal/mol) for the Finkelstein reaction (Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−}) obtained with the mean-field QM/MM method. Δ*A* ^{‡} is calculated as *A*(ξ = 0.0) − *A*(ξ = −4.0). The solvents are described with the polarizable CRK model. Values in parentheses are obtained with the nonpolarizable solvent models. The basis sets are 6-311++G(3*df*,3*p*) for MP2 and 6-31+G(*d*,*p*) for BHHLYP method.

Free energy barrier Δ*A* ^{‡} (in kcal/mol) for the Finkelstein reaction (Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−}) obtained with the mean-field QM/MM method. Δ*A* ^{‡} is calculated as *A*(ξ = 0.0) − *A*(ξ = −4.0). The solvents are described with the polarizable CRK model. Values in parentheses are obtained with the nonpolarizable solvent models. The basis sets are 6-311++G(3*df*,3*p*) for MP2 and 6-31+G(*d*,*p*) for BHHLYP method.

Free energy correction for statistical fluctuations of the QM wave function [namely, the second term in Eq. (C1)] calculated for the Finkelstein reaction in solution (in kcal/mol). The QM level is BHHLYP/6-31+G(*d*,*p*) and the solvents are described with the polarizable CRK model. Values in parentheses are obtained with the nonpolarizable solvent models.

Free energy correction for statistical fluctuations of the QM wave function [namely, the second term in Eq. (C1)] calculated for the Finkelstein reaction in solution (in kcal/mol). The QM level is BHHLYP/6-31+G(*d*,*p*) and the solvents are described with the polarizable CRK model. Values in parentheses are obtained with the nonpolarizable solvent models.

Activation free energy Δ*G* ^{‡} (in kcal/mol) for the Finkelstein reaction Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−} calculated with the mean-field QM/MM method. The QM calculation is performed at the BHHLYP/6-31+G(*d*,p) level, and the solvents are described with the polarizable CRK model. The values include the corrections for solute thermal motions (4.3 kcal/mol, see the main text) and statistical fluctuations of the QM wave function (Table II). Values in parentheses are the results obtained with the nonpolarizable models. refers to the experimental estimate.^{85}

Activation free energy Δ*G* ^{‡} (in kcal/mol) for the Finkelstein reaction Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−} calculated with the mean-field QM/MM method. The QM calculation is performed at the BHHLYP/6-31+G(*d*,p) level, and the solvents are described with the polarizable CRK model. The values include the corrections for solute thermal motions (4.3 kcal/mol, see the main text) and statistical fluctuations of the QM wave function (Table II). Values in parentheses are the results obtained with the nonpolarizable models. refers to the experimental estimate.^{85}

Dielectric constants and molecular polarizability of water, methanol (MeOH), acetonitrile (MeCN), acetone, N,N-dimethylformamide (DMF), and cyclohexane (CHX).

Dielectric constants and molecular polarizability of water, methanol (MeOH), acetonitrile (MeCN), acetone, N,N-dimethylformamide (DMF), and cyclohexane (CHX).

Density ρ (in g/cm^{3}) and vaporization enthalpy Δ*H* _{ v } (in kcal/mol) of bulk solvents calculated with the CRK and empirical (nonpolarizable) models. Vaporization enthalpy was calculated as Δ*H* _{ v } ≃ −⟨*E*⟩ + *RT*, where ⟨*E*⟩ is the average interaction energy per molecule and *R* is the gas constant.

Density ρ (in g/cm^{3}) and vaporization enthalpy Δ*H* _{ v } (in kcal/mol) of bulk solvents calculated with the CRK and empirical (nonpolarizable) models. Vaporization enthalpy was calculated as Δ*H* _{ v } ≃ −⟨*E*⟩ + *RT*, where ⟨*E*⟩ is the average interaction energy per molecule and *R* is the gas constant.

Density ρ (in g/cm^{3}) and vaporization enthalpy Δ*H* _{ v } (in kcal/mol) of bulk water and acetone at 298 K and 1 atm calculated with the polarizable and nonpolarizable models. “pol-TIP3P” refers to the CRK water model derived from the TIP3P model, while “pol-SPC” refers to the CRK water model derived from the SPC model. The CRK model for acetone is derived from the OPLS-AA model. “TIP3P,” “SPC,” and “OPLS-AA” refer to the standard empirical models. The value of Coulomb damping parameter *A* is given in parentheses. In the main text, *A* was set to 2.7 for water and 2.6 for other solvents. The Δ*A* ^{‡} stands for the free energy barrier (in kcal/mol) for the Finkelstein reaction Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−} calculated with the mean-field QM/MM method at the BHHLYP/6-31+G(*d*,*p*) level. The CRK matrix was calculated with the B3LYP/aug-cc-pVTZ method unless otherwise noted.

Density ρ (in g/cm^{3}) and vaporization enthalpy Δ*H* _{ v } (in kcal/mol) of bulk water and acetone at 298 K and 1 atm calculated with the polarizable and nonpolarizable models. “pol-TIP3P” refers to the CRK water model derived from the TIP3P model, while “pol-SPC” refers to the CRK water model derived from the SPC model. The CRK model for acetone is derived from the OPLS-AA model. “TIP3P,” “SPC,” and “OPLS-AA” refer to the standard empirical models. The value of Coulomb damping parameter *A* is given in parentheses. In the main text, *A* was set to 2.7 for water and 2.6 for other solvents. The Δ*A* ^{‡} stands for the free energy barrier (in kcal/mol) for the Finkelstein reaction Cl^{−} + CH_{3}Cl → ClCH_{3} + Cl^{−} calculated with the mean-field QM/MM method at the BHHLYP/6-31+G(*d*,*p*) level. The CRK matrix was calculated with the B3LYP/aug-cc-pVTZ method unless otherwise noted.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content