pubs.acs.org/JPCL
Explicit Solvent Model for Quantum Monte Carlo Hyung Min Cho and William A. Lester Jr* Department of Chemistry, Kenneth S. Pitzer Center for Theoretical Chemistry, University of California, Berkeley, California 94720-1460, United States, and Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, United States
ABSTRACT We present an explicit solvent model for quantum Monte Carlo (QMC) based on a hybrid quantum mechanics/molecular mechanics approach. An effective coupling Hamiltonian that combines QMC and a molecular mechanics force field is proposed to describe the interaction between QM electrons and solvent nuclei at short distances. Considering the level of approximation of the scheme, the introduction of QMC for the QM part is shown to yield significant improvement in accuracy. Results of solute-solvent interaction energy calculations on a set of hydrogen-bonded systems are found to be in excellent accord with the results of other ab initio approaches. SECTION Molecular Structure, Quantum Chemistry, General Theory
A
The first two terms account for the classical electrostatic interactions between charged particles. Here ZA is the nuclear charge of QM atom A and QR is the partial charge on solvent site R, riR is the displacement of the ith QM electron from the particle R, and RAR specifies the distance between QM atom A and solvent site R. The last term on the right-hand side of eq 2 describes the nonbonded van der Waals interactions between the QM and the MM (solvent) regions. A complication in the development of a QMC solvent model lies in the inability to use readily this typical coupling scheme. Because the QMC description involves individual electrons moving in space stochastically rather than governed by wave function electron density distributions, particular attention must be directed at developing an effective coupling Hamiltonian that modifies the electrostatic potential near the solvent atoms to eliminate nonphysical singularities and to ensure computational stability without loss of desired accuracy. In formulating the present interaction Hamiltonian, we have especially focused on the first term of eq 2 because this classical expression is not physically acceptable when a QM electron approaches a positively charged MM particle. Although this issue has been addressed in other QM/MM approaches and efforts have been made to improve the short-range behavior of the interaction potential,13-15 it leads to greater instability in the QMC method. On the basis of the examination of several coupling schemes, we propose the following interaction Hamiltonian based on a partial charge expansion15 and the exchange-repulsion
s one of the most accurate many-body electronic structure methods, quantum Monte Carlo (QMC) has been successfully applied to a wide range of molecular systems of physical and chemical interest.1-5 Compared with other first principles basis-set electronic structure methods, QMC has many attractive features. It is relatively insensitive to choice of basis set, scales naturally as N3 and with recent improvements as N2,6-9 and can accurately treat weak interactions.10-12 Furthermore, its intrinsic algorithmic adaptability to parallel computers makes the QMC method of increasing interest with advances in high-performance parallel computing technology. However, most interest today has focused on molecular systems in the gas phase, and QMC studies of chemical reactions in the condensed phase remain a largely unexplored area because of high computational cost. We present here an explicit solvent model for QMC that takes into account effectively interactions between the chemically most important active region and its aqueous surroundings. This development is achieved in a framework of a hybrid approach of quantum mechanics (QM) and molecular mechanics (MM). For a given molecular system, the effective Hamiltonian for the total system may be written, ^ solvent þ H ^ int ^ ¼ H ^ QM þ H H
ð1Þ
^ QM is the QM Hamiltonian of the active region and where H ^ solvent is a MM Hamiltonian that describes the remainder of H the system. In conventional coupling schemes, the interaction ^ int, is given by between these regions, H ^ int ¼ H
QM solvent X QR X
þ
QM solvent X ZA QR X
riR RAR R A ( ) QM solvent X X σAR 12 σ AR 6 þ 4εAR RAR RAR R A i
r 2010 American Chemical Society
R
Received Date: September 26, 2010 Accepted Date: November 9, 2010 Published on Web Date: November 15, 2010
ð2Þ
3376
DOI: 10.1021/jz101336e |J. Phys. Chem. Lett. 2010, 1, 3376–3379
pubs.acs.org/JPCL
Figure 2. Structure of selected water-solvated complexes. (a) Tetrahedral water pentamer, (b) HCONH2 þ H2O, (c) HCHO þ H2O, (d) CH3CN þ H2O, and (e) formamide solvated by eight water molecules.
Waals interaction term, we use the CHARMM atomic parameters and combination rules.16 The water dimer was chosen for initial parametrization ^ int, of and assessment of the proposed interaction potential, H eq 3. One of the water molecules is treated as a QM system by the rigorous diffusion Monte Carlo (DMC) method,1 whereas the empirical TIP3P potential17 was adopted for the other (MM) water fragment. The all-electron DMC calculations were carried out using a product function of a single determinant wave function consisting of BLYP density functional molecular orbitals with a triple-ζ polarized (TZP) basis set and a correlation function. The latter consisted of functions explicitly dependent on electron-electron and electron-nucleus distances. In subsequent calculations, we used a variant of TIP3P, flexible TIP3P, with parameters from ref 18 for maximum geometrical flexibility for the solvent water molecules. A set of water dimer configurations was selected that included the global minimum structure. CCSD(T) computations with an aug-cc-pVTZ atomic basis set were then performed to establish a reference interaction energy at each configuration, which is defined as the energy difference between the water dimer complex and two noninteracting water molecules. An optimal set of parameters was generated by minimizing the difference between the QMC/MM and the reference interaction energy at the chosen configurations. Figure 1a shows the interaction energies calculated using the QMC solvent model at the dimer configurations and several intermolecular O-O separations compared with reference energies computed by the CCSD(T) method. Full DMC calculations are also listed for comparison. The interaction energies calculated by the CCSD(T)and QMC/MM methods at different rotational configurations of the water molecules at fixed O-O distance are shown in Figure 1b along with full DMC results. The resulting QMC/MM interaction energies for water dimer are in a good agreement with those from the CCSD(T)and the full QMC calculations, representing a notable improvement over those estimated by the TIP3P potential model. A further example is provided by the small water cluster shown in Figure 2a depicting a single water molecule surrounded by four other water molecules that models the first solvation shell in liquid water. It provides a test of the transferability of parameters implemented in the QMC/MM interaction Hamiltonian. In Table 1, the interaction energy between the central QM water molecule and the surrounding four MM water molecules calculated by using the QMC solvent model is compared with MP2 (second-order MoellerPlesset perturbation theory) results obtained using selected
Figure 1. (a) Water dimer interaction energy computed by the QMC/MM method compared with those obtained from various computational methods at chosen reference configurations with different intermolecular O-O separations. (b) Interaction energy computed at different rotational configurations. The angle θ (in degrees) is formed by O-O axis of the dimer and OH bond of the QM water, representing the rotational deviation from the global minimum structure. The QMC/MM and full QMC calculations are shown with associated error bars.
fit to achieve concurrently computational stability and accuracy ^ int ¼ H
QM X
solvent X
i
R
- ξR expð - 2ξR riR Þg þ þ
QM solvent X ZA QR X A
R
RAR
þ
QR
solvent X X R
QM solvent X X A
1 expð - 2ξR riR Þ riR riR
R
akR riR k - 1 expð - bkR riR ÞÞ
k
4εAR
(
σAR RAR
12
σAR RAR
6 )
ð3Þ In the original partial charge expansion method, the ξR values are parameters that control the spread of the MM partial charges about their atomic centers. For the present purpose, however, the ξR values are taken as a set of adjustable parameters that modify the electrostatic potential experienced by a QM electron near a MM atom along with the other parameters, akR and bkR. For the van der
r 2010 American Chemical Society
3377
DOI: 10.1021/jz101336e |J. Phys. Chem. Lett. 2010, 1, 3376–3379
pubs.acs.org/JPCL
Table 1. Interaction Energy between the Central Water and the Surrounding Waters in the Tetrahedral Water Pentamer (Figure 2a)a
Table 3. Interaction Energy between a Formamide Molecule and Neighboring Eight Water Molecules (See Figure 2e)a with CP correction
with CP correction MP2/6-311(2þ,2þ)G(3df,3pd)
-24.87
-21.44
MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ
-24.15 -24.25
-20.00 -22.04
MP2/aug-cc-pVQZ
-24.05
-22.91
QMC/MM
-21.9(1)
-32.04 -33.26
B3LYP/6-311þþG(2df,2pd)
-31.52
-30.33
B3LYP/aug-cc-pVTZ
-30.83
-30.44
QMC/MM
-33.0(1)
RI-MP2 and B3LYP calculations with two basis set are presented. Counterpoise (CP)-corrected results are also given to compare with QMC/MM calculation. Energies are in kilocalories per mole.
MP2 interaction energy calculations are performed with different basis sets. For QMC/MM, the central water fragment is treated by QMC, and the other waters are represented by the TIP3P potential. Energies are in kilocalories per mole.
demonstrate the general applicability of the proposed QMC solvent model in a more realistic water-solvated environment. The initial geometry of solvated formamide was obtained by simulated annealing22 using the OPLS force field,23,24 followed by further refinement at the RI-MP2/6-311þþ(2df,2pd) level of theory. Because the present focus is to assess the performance of the present QMC water model and not to locate the global minimum structure, no further effort was made to explore configuration space to find the global minimum. The interaction energies calculated by the RI-MP2 and B3LYP methods with two basis sets are presented in Table 3 and compared with values computed using the present QMC/MM model. Even with relatively large basis sets, the RI-MP2 computed results without CP correction show ∼10% deviations from the CP-corrected counterparts. The B3LYP interaction energies are less sensitive to BSSE but show noticeable differences from the RI-MP2 and QMC/MM computed results. The QMC/MM interaction energies are in good accord with those calculated using the RI-MP2 method with CP correction. To summarize, we have proposed an explicit solvent model for QMC based on a hybrid QM/MM approach. The main obstacle in the development of the QMC/MM coupling scheme was the elimination of singularities between QM electrons and MM atoms that arise in alternative QM/MM formalisms. The modification introduced here to address this issue takes into account approximately a missing interaction in typical QM/MM implementations, namely, the exchangerepulsion experienced by individual QM electrons as they approach MM atoms. We have shown that the solute-solvent interaction energies computed using the present QMC/MM coupling scheme are consistently in good agreement with those computed by other high-level ab initio calculations, confirming good accuracy of the current scheme as a solvent model for QMC and the transferability of parameters in the coupling potential with the TIP3P water potential model for the MM part. Considering the level of approximations of the present approach, introducing QMC for the QM approach represents a major improvement in accuracy. Note also that the effect of electronic polarization in the QM treatment is intrinsically accounted for by adopting the DMC method for the QM approach, leading to greater accuracy. Further improvement is anticipated with the inclusion of solvent polarization. We are currently exploring approaches to incorporate the mutual polarization between solute and solvent regions selfconsistently by the introduction of a polarizable water potential.
Table 2. Interaction Energies of Selected Hydrogen-Bonded Complexes (Figure 2b-d)a RI-MP2/aug-cc-pVQZ
QMC/MM
-10.06
-9.6(1)
HCHO þ H2O
-5.39
-5.4(2)
CH3CN þ H2O
-4.90
-5.5(1)
a RI-MP2 energies are counterpoise (CP)-corrected. Energies are in kilocalories per mole.
basis sets. The counterpoise (CP) correction19,20 was applied to the latter to compensate for basis set superposition error (BSSE)19,21 that was found to be quite large even for the relatively large basis sets used here. Note that the present QMC/MM approach is in good agreement with CP-corrected MP2 results. The effective range of adjustment of the electrostatic interaction between QMC electrons and the MM region should be well-localized near the MM solvent nuclei so that the classical electrostatic interaction potential is recovered when QM electrons are far from MM nuclei. Because QM electrons remain essentially within the first solvation shell, the interaction between solute and solvent molecules, separated by distances larger than the range of the first solvation shell, is well represented by the classical interaction potential characterized by Coulomb interactions between the QM electrons and atomic partial charges of the MM nuclei. The number of reference systems needed to parametrize the QMC/MM coupling potential is reduced greatly owing to the relatively short distance of the effective range of adjustment from the QM solute molecule. Also, the parameters controlling this region are determined primarily by the characteristics of the MM atoms and not the QMC electrons, which implies good parameter transferability among solute systems. The performance of the proposed coupling scheme as an explicit QMC solvent model is examined further on the small representative set of hydrogen-bonded systems shown in Figure 2b-d). The computed QMC/MM interaction energies are presented in Table 2 and compared with CP-corrected RI-MP2/aug-cc-pVQZ results. As expected, the QMC/MM interaction energies yield good agreement with other ab initio methods validating parameter transferability. A larger water-solvated complex, a formamide molecule solvated by eight waters, as shown in Figure 2e, is examined to
r 2010 American Chemical Society
-36.12 -35.99
a
a
HCONH2 þ H2O
RI-MP2/6-311þþG(2df,2pd) RI-MP2/aug-cc-pVTZ
3378
DOI: 10.1021/jz101336e |J. Phys. Chem. Lett. 2010, 1, 3376–3379
pubs.acs.org/JPCL
The present explicit solvent model for QMC is an application of the QMC/MM approach to water-solvated systems. We emphasize that the present hybrid QMC/MM approach can be straightforwardly extended to other applications beyond the solvent model without further algorithmic development. With appropriate parametrization, the present QMC/MM approach can be a practical computational tool for studies of various chemical processes that arise in biological settings by providing highly accurate energy estimates for the active region of the system while maintaining an appropriate level of accuracy for interactions with surroundings.
(14)
(15)
(16)
(17)
AUTHOR INFORMATION Corresponding Author:
(18)
*To whom correspondence should be addressed. E-mail: walester@ lbl.gov.
(19)
ACKNOWLEDGMENT This work was supported by the Director,
(20)
Office of Science, of the U.S. Department of Energy under contract no. DE-AC02-05CH11231. The calculations were carried out in part at the National Energy Research Supercomputer Center (NERSC).
(21)
REFERENCES (1)
Hammond, B. L.; Lester, W. A., Jr.; Reynolds, P. J. Monte Carlo Methods in Ab Initio Quantum Chemistry; World Scientific: Singapore, 1994. (2) Anderson, J. B. Quantum Monte Carlo: Atoms, Molecules, Clusters, Liquids, and Solids, Reviews in Computational Chemistry; Wiley: New York, 1999; Vol. 13, pp 132-182. (3) L€ uchow, A.; Anderson, J. B. Monte Carlo Methods in Electronic Structures for Large Systems. Annu. Rev. Phys. Chem. 2000, 51, 501–526. (4) Foulkes, M.; Mitas, L.; Needs, R.; Rajagopal, G. Quantum Monte Carlo Simulations of Solids. Rev. Mod. Phys. 2001, 73, 33–83. (5) Lester, W. A., Jr.; Mitas, L.; Hammond, B. Quantum Monte Carlo for Atoms, Molecules and Solids. Chem. Phys. Lett. 2009, 17, 1–10. (6) Williamson, A. J.; Hood, R. Q.; Grossman, J. C. Linear-Scaling Quantum Monte Carlo Calculations. Phys. Rev. Lett. 2001, 87, 246406. (7) Manten, S.; L€ uchow, A. Linear Scaling for the Local Energy in Quantum Monte Carlo. J. Chem. Phys. 2003, 119, 1307–1312. (8) Aspuru-Guzik, A.; Salomon-Ferrer, R.; Austin, B.; Lester, W. A., Jr. A Sparse Algorithm for the Evaluation of the Local Energy in Quantum Monte Carlo. J. Comput. Chem. 2005, 26, 708–715. (9) Reboredo, F. A.; Williamson, A. J. Optimized Nonorthogonal Localized Orbitals for Linear Scaling Quantum Monte Carlo Calculations. Phys. Rev. B 2005, 71, 121105. (10) L€ uchow, A. Weak Intermolecular Interactions Calculated with Diffusion Monte Carlo. J. Chem. Phys. 2005, 123, 184106. (11) Sorella, S.; Casula, M.; Rocca, D. Weak Binding between Two Aromatic Rings: Feeling the Van Der Waals Attraction by Quantum Monte Carlo Methods. J. Chem. Phys. 2007, 127, 014105. (12) Sterpone, F.; Spanu, L.; Ferraro, L.; Sorella, S.; Guidoni, L. Dissecting the Hydrogen Bond: A Quantum Monte Carlo Approach. J. Chem. Theory Comput. 2008, 4, 1428–1434. (13) Eichinger, M.; Tavan, P.; Hutter, J.; Parrinello, M. A. Hybrid Method for Solutes in Complex Solvents: Density Functional Theory Combined with Empirical Force Fields. J. Chem. Phys. 1999, 110, 10452–10467.
r 2010 American Chemical Society
(22) (23)
(24)
3379
Laio, A.; VandeVondele, J.; Rothlisberger, U. A. Hamiltonian Electrostatic Coupling Scheme for Hybrid Car-Parrinello Molecular Dynamics Simulations. J. Chem. Phys. 2002, 116, 6941–6947. Biswas, P. K.; Gogonea, V. A. Regularized and Renormalized Electrostatic Coupling Hamiltonian for Hybrid QuantumMechanical-Molecular-Mechanical Calculations. J. Chem. Phys. 2005, 123, 164114. Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. Charmm: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187–217. Jorgensen, W. L.; Chandrasekhar, J.; D, M. J.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926–935. Schmitt, U. W.; Voth, G. A. The Computer Simulation of Proton Transport in Water. J. Chem. Phys. 1999, 111, 9361–9381. Jansen, H. B.; Ross, P. Non-Empirical Molecular Orbital Calculations on the Protonation of Carbon Monoxide. Chem. Phys. Lett. 1969, 3, 140–143. Boys, S. B.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553–566. Liu, B.; McLean, A. D. Accurate Calculation of the Attractive Interaction of Two Ground State Helium Atoms. J. Chem. Phys. 1973, 59, 4557–4558. Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220, 671–680. Jorgensen, W. L.; Tirado-Rives, J. The OPLS [Optimized Potentials for Liquid Simulations] Potential Functions for Proteins, Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657–1666. Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225–11236.
DOI: 10.1021/jz101336e |J. Phys. Chem. Lett. 2010, 1, 3376–3379
