Abstract
We present an approach to calculate the free energy profile along a condensed-phase reaction path based on high-level electronic structure methods for the reactive region. The bulk of statistical averaging is shifted toward less expensive descriptions by using a hierarchy of representations that includes molecular mechanics, density functional theory, and coupled cluster theories. As an application of this approach we study the reaction of CH Cl3 with O H- in aqueous solution.
| Original language | English |
|---|---|
| Article number | 051102 |
| Journal | Journal of Chemical Physics |
| Volume | 127 |
| Issue number | 5 |
| DOIs | |
| Publication status | Published - 2007 |
| Externally published | Yes |
ASJC Scopus subject areas
- General Physics and Astronomy
- Physical and Theoretical Chemistry
Access to Document
Other files and links
Fingerprint
Dive into the research topics of 'Hybrid approach for free energy calculations with high-level methods: Application to the SN 2 reaction of CH Cl3 and O H- in water'. Together they form a unique fingerprint.Cite this
- APA
- Standard
- Harvard
- Vancouver
- Author
- BIBTEX
- RIS
In: Journal of Chemical Physics, Vol. 127, No. 5, 051102, 2007.
Research output: Contribution to journal › Article › peer-review
}
TY - JOUR
T1 - Hybrid approach for free energy calculations with high-level methods
T2 - Application to the SN 2 reaction of CH Cl3 and O H- in water
AU - Valiev, Marat
AU - Garrett, Bruce C.
AU - Tsai, Ming Kang
AU - Kowalski, Karol
AU - Kathmann, Shawn M.
AU - Schenter, Gregory K.
AU - Dupuis, Michel
N1 - Funding Information: C H Cl 3 with O H − Valiev Marat a) Molecular Sciences Software Group, Environmental Molecular Sciences Laboratory , Pacific Northwest National Laboratory, Richland, Washington 99352, USA Garrett Bruce C. Tsai Ming-Kang Chemical and Materials Sciences Division, Pacific Northwest National Laboratory , Richland, Washington 99352, USA Kowalski Karol Molecular Sciences Software Group, Environmental Molecular Sciences Laboratory , Pacific Northwest National Laboratory, Richland, Washington 99352, USA Kathmann Shawn M. Schenter Gregory K. Dupuis Michel Chemical and Materials Sciences Division, Pacific Northwest National Laboratory , Richland, Washington 99352, USA a) Electronic mail: [email protected] 07 08 2007 127 5 051102 15 05 2007 10 07 2007 07 08 2007 2007-08-07T15:23:33 2007 American Institute of Physics 0021-9606/2007/127(5)/051102/4/ $23.00 We present an approach to calculate the free energy profile along a condensed-phase reaction path based on high-level electronic structure methods for the reactive region. The bulk of statistical averaging is shifted toward less expensive descriptions by using a hierarchy of representations that includes molecular mechanics, density functional theory, and coupled cluster theories. As an application of this approach we study the reaction of C H Cl 3 with O H − in aqueous solution. ONR DOE DOE There are many application areas that can greatly benefit from accurate electronic structure description of reaction processes provided by high-level methods, such as coupled cluster (CC) theory. 1 However, in many practical cases, explicit consideration of the complex heterogeneous environment renders such applications nearly intractable because of the need to sample over a large number of degrees of freedom. Previously we have shown 2 how to incorporate high-level ab initio methods for calculations of finite temperature averages of excited states. In this work we describe a systematic approach for CC-based calculations of a free energy profile along a reaction path—a quantity central to the analysis of reactions in large complex chemical systems. Similar in spirit to the ONIOM method, 3 our approach combines CC, density functional theory (DFT), and molecular mechanics (MM) methods to construct a sequence of multiple representations of the reactive region allowing the bulk of statistical averaging to be shifted toward less expensive descriptions. Let us consider the combined solute and solvent system. By construction, the solute contains all the reactive species requiring a quantum-mechanical (QM) description. Adopting a classical description for the solvent, the total energy of the system in the combined quantum mechanical and molecular mechanics (QM/MM) approach 4 can be expressed as a sum of quantum (QM) and classical molecular mechanics (MM) contributions, E = E QM ( r , R ; ψ ) + E MM ( r , R ) , (1) where ψ represents the electronic degrees of freedom corresponding to the ground state, and r and R are solute and solvent coordinates, respectively. The QM energy can be further decomposed into internal and external parts, E QM = E QM int ( r ; ψ ) + E QM ext ( r , R ; ρ ) . (2) Formally the internal part E QM int coincides with the solute QM energy in the gas phase. In a CC-based description, it can be represented as E CC int = ⟨ Φ ∣ H e T ∣ Φ ⟩ , (3) where H is a standard many-electron Hamiltonian and Φ is a reference many-electron wave function. The external part of the QM energy contains electrostatic interactions between solute electron density ( ρ ) and solvent classical charges ( Z I ) , E QM ext ( r , R ; ρ ) = ∑ I ∫ Z I ρ ( r ′ ) ∣ R I − r ′ ∣ d r ′ . (4) The van der Waals and Coulomb nuclear solute-solvent interactions are absorbed into the second term in Eq. (1) , which also contains the classical energy of the solvent. We assume that the reaction can be fully characterized by a reduced set of coordinates, involving only the solute part of the system. Using an equilibrium solvation model 5 the solvent degrees of freedom can be integrated out leading to the definition of the effective thermodynamic potential, W ( r , β ) = − 1 β ln ∫ e − β E ( r , R ; ψ ) d R , (5) where β = 1 ∕ k T . In many ways the effective thermodynamic potential is reminiscent of the total energy quantity, which facilitates the reaction analysis in terms of concepts developed for gas phase reaction calculations. Prior to discussing the calculation of W ( r , β ) along the reaction pathway in solute coordinate space, let us briefly address an issue of how such a pathway could be generated in the first place. Using gas phase calculations for this purpose may not always be appropriate as the solvent may significantly influence the reaction process. The main challenge is that the direct calculation of the derivatives of W ( r , β ) required for path calculations is impractical with CC or even DFT methods because of numerous electronic structure calculations during solvent averaging. The computations, however, can be significantly simplified in the zero temperature limit. In this case the presence of the solvent is taken into account, but the influence of finite temperature fluctuations on the structure of reaction path is neglected. In the zero temperature limit the major contribution to the derivatives of W ( r , β ) comes from the region where E ( r , R ; ψ ) is near its minimum with respect to R , therefore lim β → ∞ ∂ r W ( r , β ) = ∂ r E ( r , R * ; ψ ) , (6) where R * is a solvent configuration minimizing E ( r , R ; ψ ) . For a large number of solvent degrees of freedom, many optimization steps may be required to determine R * . As the solute region is kept fixed during this procedure, the computations can be simplified by following the previously suggested method 6 of introducing an effective classical representation for the solute-solvent interactions, E QM ext = ∑ i , I Z I Q i ∣ R I − r i ∣ ≡ E ESP ( r , R ; Q ) . (7) Here the QM atoms are represented by the effective electrostatic potential (ESP) charges ( Q i ) such that the electrostatic potential outside the solute region is the same as that produced from the full electron density ρ ( r ) . The charges Q i do depend on the solvent coordinates, but as long as the optimization cycle is convergent the frequency of their updates is largely irrelevant and can be exploited to achieve the best computational efficiency. The approach presented here allows straightforward incorporation of various reaction pathway methodologies. For example, the nudged elastic band (NEB) method 7 can be implemented in the following way. After generation of an initial trial pathway and determination of effective charges for each of the resulting replicas of the QM region, the self-consistent cycle is carried out consisting of (i) optimization of the solvent MM region in the field of the effective charge representation of the QM region, (ii) calculation of new forces and effective charges for the QM region through the solution of the Schrodinger equation, and (iii) evolution of the QM region according to the NEB methodology. 7 The cycle is repeated until convergence is achieved. This methodology can be applied in conjunction with both CC or less expensive DFT versions of QM/MM methods. The free energy difference between points r A and r B on the reaction path can be written as Δ W A B = − 1 β ln ⟨ e − β ( E ( r B , R ; ψ B ) − E ( r A , R ; ψ A ) ) ⟩ r A , (8) where angular brackets denote a statistical average over solvent configurations with fixed solute geometry, ⟨ … ⟩ r = ∫ … e − β E ( r , R ; ψ ) d R ∫ e − β E ( r , R ; ψ ) d R . (9) Direct QM/MM evaluation of Eq. (8) is impractical because of numerous electronic structure calculations during statistical averaging. Taking advantage of a state function property of W ( r , β ) , an alternative thermodynamic cycle can be devised where the burden of statistical sampling is shifted to a less expensive representation of the solute-solvent interactions. 6,8 For a CC/MM based description, we suggest a multitiered mapping involving three different representations of varied complexity: CC/MM, DFT/MM, and ESP/MM (see Fig. 1 ). The ESP/MM representation is obtained by replacing E QM in Eq. (1) by the effective charge energy E ESP [see Eq. (7) ]. The total free energy difference Δ W A B can then be represented as a sum of three contributions (see Fig. 1 ) Δ W A B = ( Δ W A A CC → DFT − Δ W B B CC → DFT ) + ( Δ W A A DFT → ESP − Δ W B B DFT → ESP ) + Δ W A B ESP . (10) The first and second terms represent free energy difference for changing the description of the fixed solute region from the CC to DFT representations and from DFT to a classical ESP representations [see Eq. (7) ], respectively. The third term represents the free energy difference for changing solute configuration from ( r A , Q A ) to ( r B , Q B ) within classical ESP/MM description. It is this last leg of the cycle that bears the largest burden of statistical sampling, which can now be accommodated easily. Since the position of the solute region is fixed during the change of representations (CC to DFT, DFT to ESP) it is reasonable to assume that there is a sufficient overlap between accessible solvent phase regions. In this case, free energy perturbation theory can be used to calculate corresponding free energy differences, Δ W A A CC → DFT = − 1 β ln ⟨ e − β Δ E A A CC → DFT ⟩ DFT ∕ MM , (11) Δ W A A DFT → ESP = − 1 β ln ⟨ e − β Δ E A A DFT → ESP ⟩ ESP ∕ MM , (12) where (13) Δ E A A CC → DFT = E CC ( r A , R ; ψ A ) − E DFT ( r A , R ; ψ A ) , Δ E A A DFT → ESP = E DFT ( r A , R ; ψ A ) − E ESP ( r A , R ; Q A ) . The benefit of using an intermediate DFT description is that in the ground state the DFT electron density should be fairly close 9 to the nearly exact density generated by CC. Since the electron density is the sole coupling parameter between the solute electronic degrees of freedom and the solvent [see Eq. (4) ], the derivatives δ E CC ∕ δ R and δ E DFT ∕ δ R produce only density- (not wave function) dependent terms and therefore δ Δ E A A CC → DFT δ R ≈ 0 ⇒ Δ W A A CC → DFT ≈ Δ E A A CC → DFT . (14) In the evaluation of Δ W A A DFT → ESP we can follow a resampling strategy 10 and evaluate Δ E A A DFT → ESP at a certain interval over the solvent trajectories expediently generated at the ESP/MM level of theory. There are several ways to evaluate the classical Δ W A B ESP term. Based on the smooth transformation between A and B configurations, (15) r λ = ( 1 − λ ) r A + λ r B . Q λ = ( 1 − λ ) Q A + λ Q B , Δ W A B ESP can be calculated as a sum of differences, Δ W A B ESP = − ∑ i 1 β ln ⟨ e − β Δ E λ i → λ i + 1 ESP ⟩ λ i , (16) where Δ E λ i → λ i + 1 ESP is the ESP energy difference between λ i + 1 and λ i windows that partition λ ∊ [ 0,1 ] interval. We implemented the above approach into NWCHEM (Ref. 11 ) and applied it to the S N 2 reaction of C H Cl 3 and O H − in aqueous solution. Only gas phase studies of this reaction have been performed 3,12 so far, and there is some uncertainty whether the S N 2 mechanism is a viable option for base-catalyzed hydrolysis of C H Cl 3 in aqueous solution. Our solute QM region contained C H Cl 3 and O H − molecules. At the DFT level it was described using the B3LYP exchange-correlation functional with the 6 - 31 + G * basis set, and the CC description was based on the locally renormalized variant of CCSD(T)(LR-CCSD(T)) (Ref. 13 ) with the aug-cc-pVDZ (Ref. 14 ) basis set. The van der Waals parameters for the QM region were taken from standard Amber force field definitions. 15 The solute was embedded into a 30 Å cubic box of 888 SPC/E (Ref. 16 ) solvent water molecules representing the MM region. Reaction pathway calculations were performed using the DFT/MM representation of the system using the zero temperature methodology described earlier. First we determined the transition state structure, characterized by a single imaginary mode of 238 cm − 1 . Then we optimized reactant and product states with the initial structures generated through the displacement of the transition state in the direction of the negative mode. To ensure full relaxation of the solvent, optimization of reactant, transition, and product states were performed in combination with molecular dynamics equilibration ( 40 ps at T = 298.15 K ) of the solvent region. The initial pathway for NEB calculations was calculated by linear interpolation between the reactant, transition, and product states of the entire solute-solvent system with six beads/replicas for each segment. In the reactant state the O H − group is positioned 3.1 Å away from the carbon on the C H Cl 3 group serving as a hydrogen bond acceptor to chloroform. During the reaction it makes a symmetrical approach to chloroform such that the OH bond bisects one of the Cl–C–Cl angles. The water solvation shell of O H − changes from six-fold coordination in the reactant state to five-fold and three-fold coordination in the transition and product states, respectively. In the product state the Cl − ion, located 3.23 Å from the C H Cl 3 O H complex, is coordinated by four waters. 17 The free energy profile, including individual contributions, is shown in Fig. 1 . For Δ W A B ESP calculations we used effective charges generated during the NEB procedure. The λ -summation interval was partitioned in 0.1 Å increments ( λ i = 0 , 0.1 , … , 1.0 ) . Using a double wide sampling strategy, 18 the free energy differences for λ i → λ i − 1 and λ i → λ i + 1 were obtained simultaneously by sampling at λ i over 20 ps of molecular dynamics simulation at constant temperature ( T = 298.15 K ) . The calculation of Δ W DFT → ESP was done by resampling the 100 ps solvent trajectory generated at the ESP/MM level with the interval of 0.5 ps . The free energy difference Δ W CC → DFT was calculated according to Eq. (14) over the NEB reaction path. To estimate the accuracy of this approximation we also performed a series of CC and DFT calculations over the 10 ps solvent trajectory with 0.5 ps interval for the reactant and transition states. The calculated standard deviations for Δ E A A CC → DFT ( 0.45 kcal ∕ mol for reactant and 0.13 kcal ∕ mol for transition states) and the estimated barrier error of 0.5 kcal ∕ mol support the validity of our approximation. Overall we find that the CC/MM approach gives a free energy reaction barrier of 29.3 kcal ∕ mol and an overall reaction free energy of − 46.7 kcal ∕ mol . We note that DFT/MM treatment of the same process underestimates the reaction barrier ( 24.6 kcal ∕ mol ) but agrees well on the reaction free energy ( − 47 kcal ∕ mol ) . In conclusion we have developed an efficient and systematic approach for free energy calculations along a condensed-phase reaction path based on high-level electronic structure methods for the reactive region. The approach can be used in conjunction with not only QM/MM but also ONIOM (Ref. 3 ) and other methodologies. 19 We have illustrated our approach by studying the S N 2 reaction of C H Cl 3 with O H − in aqueous phase. Significant differences between free energy reaction barriers between CC and DFT descriptions illustrate the importance of an accurate description of electron correlation effects in free energy calculations. In agreement with experiment, 20 the high reaction barrier obtained in our calculation indicates that the S N 2 mechanism is unlikely to play a major role in the hydrolysis of C H Cl 3 in aqueous solution. Three of the authors (M.V., B.C.G., and S.M.K.) would like to acknowledge support from Office of Naval Research. Three of the authors (M.K.T., G.K.S., and M.D.) would like to acknowledge support from the Division of Chemical Sciences, Office of Basic Energy, DOE. This research was performed in part using the MSCF in EMSL, a national scientific user facility sponsored by the U.S. DOE, OBER and located at PNNL. Battelle operates PNNL for DOE. FIG. 1. (Color online) Thermodynamic cycle for free energy calculations (top left); free energy profile (bottom left) and reactant, transition, and product states (right) for the reaction C H Cl 3 + O H − → C H Cl 2 O H + Cl − .
PY - 2007
Y1 - 2007
N2 - We present an approach to calculate the free energy profile along a condensed-phase reaction path based on high-level electronic structure methods for the reactive region. The bulk of statistical averaging is shifted toward less expensive descriptions by using a hierarchy of representations that includes molecular mechanics, density functional theory, and coupled cluster theories. As an application of this approach we study the reaction of CH Cl3 with O H- in aqueous solution.
AB - We present an approach to calculate the free energy profile along a condensed-phase reaction path based on high-level electronic structure methods for the reactive region. The bulk of statistical averaging is shifted toward less expensive descriptions by using a hierarchy of representations that includes molecular mechanics, density functional theory, and coupled cluster theories. As an application of this approach we study the reaction of CH Cl3 with O H- in aqueous solution.
UR - http://www.scopus.com/inward/record.url?scp=34547685930&partnerID=8YFLogxK
UR - http://www.scopus.com/inward/citedby.url?scp=34547685930&partnerID=8YFLogxK
U2 - 10.1063/1.2768343
DO - 10.1063/1.2768343
M3 - Article
AN - SCOPUS:34547685930
SN - 0021-9606
VL - 127
JO - Journal of Chemical Physics
JF - Journal of Chemical Physics
IS - 5
M1 - 051102
ER -