Molecular Mechanics Liquid Water

pubs.acs.org/JPCL

Compatibility of Quantum Chemical Methods and Empirical (MM) Water Models in Quantum Mechanics/ Molecular Mechanics Liquid Water Simulations Katherine E. Shaw, Christopher J. Woods, and Adrian J. Mulholland* Centre for Computational Chemistry, School of Chemistry, University of Bristol, Bristol BS8 1TS, U.K.

ABSTRACT High-level (MP2, HF, and BLYP with the aug-cc-pVDZ basis set) quantum mechanics/molecular mechanics (QM/MM) Monte Carlo free energy simulations of liquid water are used here to test the compatibility of various QM methods with four standard empirical “molecular mechanics” (MM) water models. Consistency of QM methods with water models is of particular importance, given the aqueous environment of many of the systems of interest for QM/MM modeling (e.g., biological systems). The results show that treating a single water molecule using a QM method in bulk TIP3P can induce solvent structuring consistent with experiment. The results also show that the TIP4P model is the most suitable water model of those tested for such QM/MM simulations, while the TIP5P model is not well suited. The findings have important implications for future QM/MM method development and applications. They indicate that the choice of MM models should be made carefully for consistency and compatibility in QM/MM simulations. SECTION Molecular Structure, Quantum Chemistry, General Theory

here are many empirical, atomistic “molecular mechanics” (MM)-type models designed to represent liquid water in molecular simulations.1 These have typically been developed to try to reproduce, as closely as possible, properties of bulk liquid water, with many of these models giving good results under specific conditions.2 These models have also provided molecular-level insight into the unique properties of liquid water. One popular method used in biomolecular simulations3 is the TIP3P model.4,5 This simple model treats H2O as a rigid three-point system, with point charges on each of the atoms and Lennard-Jones parameters only on the oxygen. However, TIP3P does not accurately reproduce the second solvation shell seen experimentally in liquid water: it shows very little, if any, structure beyond the first solvation shell.2,6,7 There is also a modified version, TIPS3P,8 used in CHARMM (e.g., CHARMM22 and CHARMM27 forcefields),9,10 which has further Lennard-Jones parameters on the hydrogen atoms. The more sophisticated TIP4P5 and TIP5P6 models use more interaction sites (4 or 5, respectively) and show a clear second solvation shell.2,5 Combined quantum mechanics/molecular mechanics (QM/MM) methods are increasingly widely used, for example, for simulations of chemical reactions in condensed phases.11-13 It is now possible to carry out QM/MM calculations with highly accurate correlated ab initio QM methods.14 Increasingly, QM/MM approaches are being applied in molecular simulations (Monte Carlo (MC) or molecular dynamics) of condensed phases, e.g., of biomolecular systems.3 Many QM/ MM calculations use “off the shelf”, standard MM parameters. This has proved reasonable in many applications,11,12 but

T

r 2009 American Chemical Society

there are clearly limitations and potential pitfalls in this ad hoc approach. We have recently developed a method15 that allows us to use free energy calculations to investigate the compatibility of QM and MM methods. In the work presented here, thermodynamic integration free energy simulations over a lambda (λ) coordinate are used as a rigorous test of QM/MM methods. The compatibility of the four different standard MM water models mentioned above is tested with three levels of QM theory: density functional theory with the BLYP functional; and ab initio Hartree-Fock (HF) and second-order Møller-Plesset perturbation theory (MP2). These methods, while not an exhaustive set, provide a good variety of QM treatments, spanning the range of QM methods typically used in current high level QM/MM modeling. Consistency of QM methods with water models is of particular importance, given the aqueous environment of many of the systems of interest for QM/MM modeling (e.g., biological systems). It is by no means guaranteed a priori that water models are suitable for QM/MM simulations. A QM method and MM model may be defined as being compatible if they are interchangeable; therefore it should be possible to replace a MM water molecule in bulk MM solvent with a QM molecule without significant structural or energetic changes. Radial distribution functions (RDFs) can be used to examine the structural effects of replacing a single MM water molecule with a QM water molecule, in bulk MM water. Received Date: October 6, 2009 Accepted Date: November 9, 2009 Published on Web Date: November 17, 2009

219

DOI: 10.1021/jz900096p |J. Phys. Chem. Lett. 2010, 1, 219–223

pubs.acs.org/JPCL

Figure 1. RDFs for different lambda (λ) values for the perturbation of a QM water (λ = 0), treated by MP2/aug-cc-pVDZ, into a MM water molecule (λ = 1) in bulk TIP3P, TIPS3P, TIP4P, and TIP5P solvent. Each λ value plot is shifted up by 0.25 to make the distinctions between the curves more obvious. Table 1. Free Energies (in kcal 3 mol-1) for the Perturbation of a QM Water Molecule into a MM Water Molecule in Bulk MM Solvent (All QM Calculations Used the aug-cc-pVDZ Basis Set)a

The oxygen-oxygen RDF between the water molecule being perturbed from QM to MM, and all surrounding bulk MM water molecules, was calculated and averaged for each value of λ. Figure 1 shows the RDFs for each of the different MM water models as a function of λ for the MP2 QMfMM perturbation. It is immediately evident that the TIP5P MP2/ MM simulation is unstable: the RDF shows that the bulk water is highly disordered for lambda values less than 0.429 (i.e., up to about 50% QM character). Similar instability has been observed before.15 Examination of the structures here shows that the lone pairs on the TIP5P molecules become embedded in the electron cloud of the QM water, a clearly unphysical geometry. The hydrogen bonding network around the QM molecule is disrupted, as can be seen in the RDF: for example the first peak for lambda values up to 0.429, is at a shorter distance. There is no detectible long-range reorganization of the bulk solvent during the perturbation. Similar problems were observed with all three QM methods. It is important to stress that this does not indicate any weakness or problem with the TIP5P model, which performs very well for MM simulations of liquid water.6 It should be remembered that none of the water models tested here were developed for QM/ MM simulations. Some water models may be more suitable than others for use in QM/MM approaches. In the TIP3P and TIPS3P simulations, a clear second solvation shell emerges as the QM character of the perturbed water increases: it is apparent even with only approximately 50% QM character (up to λ = 0.429). The peak of the first solvation shell falls at approximately 2.75 Å for all λ values, which is consistent with experiment (X-ray scattering and neutron diffraction experiments under ambient conditions place the peak between 2.73 and 2.88 Å16-19). It is interesting that the introduction of a single QM water molecule in bulk TIP3P solvent induces this important structural feature (i.e., the second solvation shell). Only the TIP4P simulations


perturbation

QM method

QM f TIP3P

BLYP

0.5 ((0.3)

HF

3.3 ((0.5)

QM f TIPS3P

QM f TIP4P

MP2

1.7 ((0.4)

BLYP

0.5 ((0.3)

HF MP2

2.7 ((0.4) 1.5 ((0.5)

BLYP

-0.1 ((0.1)

HF QM f TIP5P

a

free energy change

2.4 ((0.2)

MP2

1.2 ((0.2)

BLYP

10.4 ((1.9)

HF

12.7 ((1.8)

MP2

11.8 ((1.9)

One standard deviation is shown in brackets.

show a near constant RDF over all lambda values, with the second solvation shell maintained throughout. These structural results show that, of the water models tested here, TIP4P is the most compatible with QM/MM simulations of this type. The RDFs for the other simulations (i.e., for the perturbation from HF or BLYP QM) can be found in the Supporting Information. They exhibit very similar trends and so lead to the same conclusions, i.e., TIP4P is also the best water model for QM/MM simulations at the BLYP and HF levels. The free energy changes for perturbing the QM water molecule into a MM water are shown in Table 1. These results also give an indication of the compatibility of each QM method with each water model. If a QM and a MM method are compatible, the free energy change would be expected ideally to be near zero. It is clear from Table 1 that the free energy

220


pubs.acs.org/JPCL

Figure 2. PMFs for the perturbation of a QM water molecule (λ = 0) into MM water (λ = 1) into TIP3P, TIPS3P, TIP4P, and TIP5P solvent. In all cases, the aug-cc-pVDZ basis set was used. The results with the TIP5P water model are shown separately because the free energies are significantly larger. The other three plots, grouped by QM method used, are all plotted on the same scale to facilitate comparison of the free energy changes and curvature.

excellent compatibility. Of the QM methods, BLYP appears to be the most compatible with TIP4P, and HF the least (the HF to TIP4P free energy change is much larger (2.4 kcal 3 mol-1)). Excessive concentration of charge in the HF calculations (which lack electron correlation) may exacerbate any problems of QM/MM overpolarization. The potentials of mean force (PMFs), generated by plotting the free energy change as a function of lambda, for these perturbations are shown in Figure 2. It is apparent from inspection that these PMFs show different levels of curvature. For example, the TIP3P and TIPS3P PMFs are more curved than for TIP4P. This difference in curvature is most likely due to the structural changes seen during the TIP3P calculation, as evidenced in the RDFs. This may be due to the (probably entropy dominated) cost of reorganizing the solvent to form a second solvation shell. The results here show that care must be taken when combining MM and QM models in QM/MM simulations. Some standard, well-parametrized MM models (such as TIP5P) may be unsuitable for QM/MM simulations. The approach applied here provides a route to testing the compatibility of QM and MM methods. The RDFs and the free energy changes (e.g., here for the perturbation of a QM water molecule into a MM water) give a good indication of the structural and energetic compatibility of each MM model with different QM methods. The lower free energy changes seen in the TIP4P simulations in conjunction with the relatively small changes in RDFs during the QM to MM perturbation are evidence that TIP4P is significantly more compatible with the current QM/MM approaches than the other (MM) water models studied here. This is true at different several levels of QM theory, using either ab initio or density functional theory (DFT)-based QM methods. The BLYP results indicate that this QM method in particular shows excellent compatibility with TIP4P.

changes, as water is perturbed from QM to MM, are all positive, except for the BLYP to TIP4P simulation. The TIP5P perturbations exhibit large free energy changes for all the QM methods (greater than 10 kcal 3 mol-1, Table 1). This is significantly larger than for the other MM water models (for which all the free energy changes are between 0 and 4 kcal 3 mol-1, Table 1). The high free energy changes found with TIP5P are likely to be the result of the severe incompatibility between TIP5P and all the QM methods, shown in the RDFs. The TIP3P and TIPS3P perturbations also display higher free energy changes than the TIP4P simulations. Part of the free energy change for TIP3P and TIPS3P can be attributed to the structural reorganization seen during the simulation. This reorganization “cost” may be largely entropic, although this has not been investigated. The entropic component of the solvent reorganization energy is likely to be larger for TIP3P than TIP4P (based on the increased ordering seen in the RDFs, Figure 1), causing this larger free energy change. The TIP4P model is clearly the most compatible of the water models tested here with all the QM methods, as demonstrated by the smaller free energy changes and less change in solvent structure. Different strengths of water-water interaction energies between the QM method and the MM water model will also play a part in determining the free energy change. The free energy changes from QM to MM are positive in all but one case (for the BLYP to TIP4P perturbation, which gives an energy close to zero); this indicates that, in almost all cases, the QM water molecule is more strongly solvated than its MM equivalent. A factor that should be considered is possible overpolarization of the QM water molecule by its MM environment (the MM models are parametrized to include polarization implicitly, while the QM method models electronic polarization explicitly).20 The free energy change for the BLYP to TIP4P perturbation is very small (-0.1 kcal 3 mol-1), indicating


221


pubs.acs.org/JPCL

Most current biomolecular QM/MM simulations apply the TIP3P or TIPS3P water model, but this may not be the best choice. Of course, in an inhomogeneous system (e.g., a protein in solution), where a variety of groups are treated by MM, consistency of the water model with the MM force field is also an essential factor. It is vital that any modification of MM parameters to improve QM/MM models should not significantly worsen the treatment of MM-MM interactions. One focus of future development work should be to test possible modifications of MM models (for water, and other molecules) to investigate and improve their compatibility with QM methods for QM/MM simulations. This could be by modification of MM (e.g., water) models (which could include changes in atomic charges or Lennard-Jones parameters, more sophisticated electrostatic descriptions beyond atomic point charge models, and MM treatment of electronic polarization) and/or QM/MM interaction schemes. Within the current QM/MM implementation, an obvious improvement is optimizing the MM Lennard-Jones parameters used for the QM atoms. Previous optimizations of this type have resulted in significant changes only for parameters on polar hydrogens.21-27 It is worth noting that, for the water models here, except for TIPS3P, there are only Lennard-Jones parameters on the oxygen. Therefore, to improve the compatibility of these models, it may be necessary to include extra parameters on more of the interaction sites in the system. For example, for TIP5P it may be necessary to include parameters on the lone pairs to prevent them from becoming embedded in the hydrogens of the QM molecule. Care would need to be taken that any reparameterizations do not adversely affect the simulated properties of bulk liquid water. Now that it is possible to perform QM/MM calculations, including molecular simulations,15 with high levels of QM theory,14 an important aspect of QM/MM method development and testing should focus on the compatibility of QM and MM models and potentially the development of more sophisticated and compatible MM models (e.g., force fields for proteins and other biomolecules). Testing of MM water models is an essential first stage in this process. Free energy simulations, of the type outlined here, will help in these developments. Ideally, each exact combination of QM and MM methods should be tested, although DFT-based QM/MM methods are generally not very sensitive to the basis set,28 so QM/MM models for, e.g., B3LYP or BLYP may be transferable to use with different basis sets.26 The free energies presented in this paper were calculated using a method outlined previously.15 This approach combines the idea of sampling a high-level ensemble (in this case a system treated using QM/MM) using a lower level reference system (in this case using pure MM), using the MetropolisHastings29 MC algorithm. The TIP3P, TIPS3P, and TIP4P systems studied here consisted of 1680 water molecules, and the TIP5P system had 1679 molecules (resulting in box dimensions of ∼36 36 36 Å), giving the correct density for the simulation temperature (298.15 K). The systems were simulated using MC methods under periodic boundary conditions with the NPT ensemble and 15 Å molecule-based electrostatic and


nonbonded cutoffs. A single QM water molecule was perturbed into a MM water molecule in the corresponding bulk MM solvent over a λ coordinate. All molecules were treated as rigid bodies, and were allowed to rotate and translate randomly in accordance with a MC acceptance test. The gas-phase QM intramolecular energy is constant, due to the molecule being treated as rigid (at the MM geometry), and is subtracted from the final QM/MM energies. The choice of geometry should have little effect on the calculated energetics. Only the intermolecular energies are sampled here. The repulsion-dispersion terms between the QM molecule and the bulk MM molecules are treated using empirical Lennard-Jones parameters. The full simulation set up is as described previously.15 In this work, eight λ values were used: 0.0, 0.142, 0.285, 0.429, 0.571, 0.714, 0.857, and 1. The energy of the systems at different λ values (corresponding to “pure” QM/MM at λ = 0, “pure” MM at λ = 1, and a hybrid of the two at intermediate values) was calculated using the equation ETotal ¼ ð1 -λÞEQM=MM þ λEMM The free energy change was calculated using replica exchange thermodynamic integration (RETI).30 In total, 10 million MC moves, per λ value, were simulated, divided into 200 iteration cycles, each consisting of 50 000 MC steps. The first 250 000 steps were not included in the analysis, to allow the system to equilibrate. All the QM calculations were run with the aug-cc-pVDZ basis set.

SUPPORTING INFORMATION AVAILABLE RDFs for the HF and BLYP QM/MM simulations. This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION Corresponding Author: *To whom correspondence should be addressed. E-mail: Adrian. [email protected].

ACKNOWLEDGMENT This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol (www.acrc.bris.ac.uk). We thank EPSRC (A.J.M. and C.J.W., Grants EP/G007705/1 and EP/G042853/1) and BBSRC (A.J.M. and K.E.S.) for support. A.J.M. is an EPSRC Leadership Fellow.

REFERENCES (1)

(2)

(3)

(4)

222

Guillot, B. A Reappraisal of What We Have Learnt during Three Decades of Computer Simulations on Water. J. Mol. Liq. 2002, 101, 219–260. Jorgensen, W. L.; Tirado-Rives, J. Potential Energy Functions for Atomic-Level Simulations of Water and Organic and Biomolecular Systems. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6665–6670. van der Kamp, M. W.; Shaw, K. E.; Woods, C. J.; Mulholland, A. J. Biomolecular Simulation and Modelling: Status, Progress and Prospects. J. R. Soc. Interface 2008, 5, S173–S190. Jorgensen, W. L. Quantum and Statistical Mechanical Studies of Liquids 0.10. Transferable Intermolecular Potential Functions


pubs.acs.org/JPCL

(5)

(6)

(7)

(8)

(9)

(10)

(11) (12) (13)

(14)

(15)

(16)

(17) (18)

(19)

(20)

(21)

for Water, Alcohols, and Ethers - Application to Liquid Water. J. Am. Chem. Soc. 1981, 103, 335–340. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926– 935. Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910–8922. Mantz, Y. A.; Chen, B.; Martyna, G. J. Structural Correlations and Motifs in Liquid Water at Selected Temperatures: Ab Initio and Empirical Model Predictions. J. Phys. Chem. B 2006, 110, 3540–3554. Neria, E.; Fischer, S.; Karplus, M. Simulation of Activation Free Energies in Molecular Systems. J. Chem. Phys. 1996, 105, 1902–1921. MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586–3616. 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. Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198–1229. Senn, H. M.; Thiel, W. QM/MM Studies of Enzymes. Curr. Opin. Chem. Biol. 2007, 11, 182–187. Lin, H.; Truhlar, D. G. QM/MM: What Have We Learned, Where Are We, and Where Do We Go from Here?. Theor. Chem. Acc. 2007, 117, 185-199. Claeyssens, F.; Harvey, J. N.; Manby, F. R.; Mata, R. A.; Mulholland, A. J.; Ranaghan, K. E.; Schutz, M.; Thiel, S.; Thiel, W.; Werner, H. J. High-Accuracy Computation of Reaction Barriers in Enzymes. Angew. Chem., Int. Ed. 2006, 45, 6856– 6859. Woods, C. J.; Manby, F. R.; Mulholland, A. J. An Efficient Method for the Calculation of Quantum Mechanics/Molecular Mechanics Free Energies. J. Chem. Phys. 2008, 128, 014109– 014101. Hura, G.; Sorenson, J. M.; Glaeser, R. M.; Head-Gordon, T. A High-Quality X-ray Scattering Experiment on Liquid Water at Ambient Conditions. J. Chem. Phys. 2000, 113, 9140–9148. Soper, A. K.; Phillips, M. G. A New Determination of the Structure of Water at 25 °C. Chem. Phys. 1986, 107, 47–60. Soper, A. K.; Bruni, F.; Ricci, M. A. Site-Site Pair Correlation Functions of Water from 25 to 400 °C: Revised Analysis of New and Old Diffraction Data. J. Chem. Phys. 1997, 106, 247– 254. Narten, A. H.; Levy, H. A. Liquid Water-Molecular Correlation Functions from X-ray Diffraction. J. Chem. Phys. 1971, 55, 2263. Senthilkumar, K.; Mujika, J. I.; Ranaghan, K. E.; Manby, F. R.; Mulholland, A. J.; Harvey, J. N. Analysis of Polarization in QM/ MM Modelling of Biologically Relevant Hydrogen Bonds. J. R. Soc. Interface 2008, 5, S207–S216. Bash, P. A.; Ho, L. L.; MacKerell, A. D.; Levine, D.; Hallstrom, P. Progress toward Chemical Accuracy in the Computer Simulation


(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

223

of Condensed Phase Reactions. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 3698–3703. Freindorf, M.; Gao, J. L. Optimization of the Lennard-Jones Parameters for a Combined Ab Initio Quantum Mechanical and Molecular Mechanical Potential Using the 3-21G Basis Set. J. Comput. Chem. 1996, 17, 386–395. Freindorf, M.; Shao, Y. H.; Furlani, T. R.; Kong, J. LennardJones Parameters for the Combined QM/MM Method Using the B3LYP/6-31þG*/AMBER Potential. J. Comput. Chem. 2005, 26, 1270–1278. Gao, J. L.; Xia, X. F. A Priori Evaluation of Aqueous Polarization Effects through Monte-Carlo QM-MM Simulations. Science 1992, 258, 631–635. Luque, F. J.; Reuter, N.; Cartier, A.; Ruiz-Lopez, M. F. Calibration of the Quantum/Classical Hamiltonian in Semiempirical QM/MM AM1 and PM3 Methods. J. Phys. Chem. A 2000, 104, 10923–10931. Pentik€ ainen, U.; Shaw, K. E.; Senthilkumar, K.; Woods, C. J.; Mulholland, A. J. Lennard-Jones Parameters for B3LYP/ CHARMM27 QM/MM Modeling of Nucleic Acid Bases. J. Chem. Theory Comput. 2009, 5, 396–410. Riccardi, D.; Li, G. H.; Cui, Q. Importance of van der Waals Interactions in QM/MM Simulations. J. Phys. Chem. B 2004, 108, 6467–6478. Hermann, J. C.; Pradon, J.; Harvey, J. N.; Mulholland, A. J. High Level QM/MM Modeling of the Formation of the Tetrahedral Intermediate in the Acylation of Wild Type and K73A Mutant TEM-1 Class A Beta Lactamase. J. Phys. Chem. A 2009, 113, 11984–11994. Hastings, W. K. Monte-Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika 1970, 57, 97–109. Woods, C. J.; Essex, J. W.; King, M. A. The Development of Replica-Exchange-Based Free-Energy Methods. J. Phys. Chem. B 2003, 107, 13703–13710.


Molecular Mechanics Liquid Water

Recommend Documents