J. Phys. Chem. B 2008, 112, 6233-6240
6233
The Polarizability of Point-Polarizable Water Models: Density Functional Theory/ Molecular Mechanics Results† Bernhard Schropp and Paul Tavan* Lehrstuhl fu¨r Biomolekulare Optik, Ludwig-Maximilians-UniVersita¨t, Oettingenstrasse 67, 80538 Mu¨nchen, Germany ReceiVed: July 21, 2007; In Final Form: October 17, 2007
Molecular dynamics (MD) simulations of bulk liquid water at different thermodynamic conditions or of biomolecules in aqueous solution require a molecular mechanics (MM) force field that accounts for the sizable electronic polarizability R of the water molecule. A considerable number of such polarizable water models has been suggested in the past. Most of them agree that one should employ the experimental value Rexp for the electronic polarizability and compute the induced dipole moment µi through linear response from the electric field E at the position rO of the oxygen atom. However, several more recent models have suggested somewhat smaller values for R. Using a hybrid method that combines density functional theory for a selected water molecule with an MM description of its liquid water environment, here we show that the choice of Rexp is solely correct if the induced dipole moment µi is calculated from the average electric field 〈E〉 within the volume occupied by the given water molecule. Because of considerable field inhomogeneities caused by the structured aqueous environment, the average field 〈E〉 is much smaller than the local spot check E(rO). However, as opposed to E(rO), the average field 〈E〉 cannot be easily calculated in MM-MD simulations. Therefore, in polarizable MM water models, one should calculate the induced dipole moment µi from E(rO) through the reduced polarizability Reff ) 0.68Rexp, which then effectively accounts for the inhomogeneities of the electric field within the volume of a water molecule embedded in liquid water.
1. Introduction Liquid water is undoubtedly the most important fluid on earth, because life originates from water. It has long been the subject of intense theoretical and experimental studies.1-3 Nevertheless, the behavior of this complex liquid still poses many puzzles, mainly because the structural and electrostatic properties of the individual water molecules strongly change as soon as they are transferred from the gas phase into the aqueous environment of the liquid phase and because these changes are not yet sufficiently well characterized.4,5 As long as water molecules are isolated in the gas phase, their structural and electrostatic properties are precisely known. Here, the geometry of a water molecule is given by the O-H bond length rOH ) 0.957 Å and the H-O-H bond angle φHOH ) 104.5°.6 The molecule carries a large electrostatic dipole moment µ ) 1.85 D.7 It has a sizable electronic polarizability,8 which, in the coordinate system displayed by Figure 1, is given by the tensor components Rxx ) 1.47 Å3, Ryy ) 1.53Å3, and Rzz ) 1.42 Å3. Thus, the electronic polarizability of the water molecule is nearly isotropic. Mainly because of this polarizability, the average molecular dipole moment is 2.1 D in the water dimer,9 increases in larger water clusters,9 and reaches, in the condensed phase, still larger values between 2.4 and 2.9 D. Here, classical molecular dynamics (MD) simulations of the dielectric properties predict values in the range 2.4-2.6 D,10,11 whereas values around 2.9 D are obtained from ab initio MD simulations12,13 and from analyses of experimental data.14,15 †
Part of the “Attila Szabo Festschrift”. * Corresponding author. E-mail:
[email protected]; phone: +49-89-2180-9220.
Figure 1. Geometry of the water molecule and coordinate system employed for its description.
The properties of isolated water molecules can be quite accurately calculated by quantum mechanics using, e.g., density functional theory (DFT).16,17 As follows from the electron density obtained by such calculations, a water molecule is nearly spherical and very small, being approximately the size of an oxygen atom. Theoretical descriptions become much more challenging if one wants to compute the properties of the bulk liquid or even of biomolecular mixtures. Because of the long-range electrostatic forces, physically realistic models of such complex systems have to be chosen very large, i.e., have to comprise at least several thousand atoms.18 The computational effort connected with such large system sizes prevents the application of precise quantum mechanical methods to the calculation of the intra- and intermolecular forces. Therefore, one has to resort to simplified molecular mechanics (MM) force fields.
10.1021/jp0757356 CCC: $40.75 © 2008 American Chemical Society Published on Web 01/17/2008
6234 J. Phys. Chem. B, Vol. 112, No. 19, 2008 MM force fields describe the water molecule by model potentials accounting for its geometry, its electrostatic signature and response, and its van der Waals interactions. The geometry is usually considered to be fixed at the experimental gas-phase finding, because at ambient temperatures the intramolecular motions are quantum mechanically frozen into the vibrational ground state and because geometry changes upon transfer into the liquid phase are expected to be small. Furthermore, the van der Waals interactions are generally represented by a single Lennard-Jones potential centered at the position of the oxygen atom. Thus, the key difference between the plethora of model potentials suggested so far for water is the treatment of the electrostatics, and each of these models represents a specific compromise between accuracy and efficiency (see the excellent review by Guillot in ref 4). The properties of the condensed phase are then calculated by MD or Monte Carlo simulations of large periodic boxes containing many thousands of MM water models and, in the case of bio-molecular mixtures, additional MM models of other molecules. As is common in MM-MD, the computational effort increases with the complexity of the respective model potential. The complexity is measured by the number of position vectors (the so-called “points”) and electrostatic moments (point charges, dipoles, etc.), which have to be handled in the computation of the interactions. Minimal water models such as the “three-point transferable intermolecular potential” (TIP3P) model of Jorgensen et al.19 or the “effective simple point charge” (SPC/E) model of Berendsen et al.20 employ only three points, which cover the positions of the three atoms, and assign three fixed partial charges to these positions. These minimal models try to account for the condensed phase effects of electronic polarization in a mean field sense by assigning partial charges in such a way that the dipole moment has an increased value of about 2.35 D. In simulations of biomolecular systems they still represent the standard. For an improved description of the higher multipole moments, which steer the local arrangements of the water molecules in the condensed phase, more complex models such as the “four- and five-point transferable intermolecular potentials” (TIP4P, TIP5P, respectively) have been suggested.19,21 These models place additional fixed partial charges to massless sites within the molecule. All these mean field models perform reasonably well on properties of pure liquid water at certain thermodynamic conditions. However, the quality of results rapidly decreases at other conditions or for inhomogeneous systems like proteinwater mixtures. Here, the water molecules are exposed to different and locally varying mean electric fields. Because of the polarizability of the molecules, any fixed dipole moment will then represent a poor approximation. By investigating a variety of molecular environments, Lamoureux et al.22 concluded that the explicit inclusion of the “polarizability is essential to get from the same water model accurate energetics in the vicinity of highly polar moieties (such as carbonyl groups), small ions (such as sodium or chloride), as well as in anisotropic nonpolar environments.” Therefore, a still increasing number of polarizable water models have been developed during the past decades. Table 1 lists a selection of polarizable water models available in the literature. The underlying concepts are summarized in a recent review by Yu and van Gunsteren.41 We have sorted these models into three main conceptual classes. The first class covers the so-called “charge-on-spring (COS)” models. Here, a large charge qc is fixed, e.g., at the position rO of the oxygen atom, and an opposite charge -qc is bound to this position through a harmonic potential with spring constant kc. The elongation of
Schropp and Tavan TABLE 1: A Selection of Polarizable Water Models from the Literaturea model
µ0x /D Rxx/Å3 Ryy/Å3 RzzÅ3
class
year 〈µx〉/D
experimenta
1.85 1.468
1.528 1.415
STR/1b COS/B1c COS/B2c COS/G2d COS/G3d SWM4-DPe SWM4-NDPf
1.95 1.90 2.07 1.85 1.85 1.85 1.85
1.445 1.401 0.930 1.255 1.250 1.043 0.978
1.445 1.401 0.930 1.255 1.250 1.043 0.978
1.445 COS 1.401 0.930 1.255 1.250 1.043 0.978
1990 2003 2003 2004 2004 2003 2006
2.73 2.65 2.57 2.53 2.52 2.41 2.37
KJg BSVh SCPDPi POL3j GCPMk
1.85 1.85 1.85 1.85 1.86
1.470 1.440 1.444 1.444 1.444
1.470 1.440 1.444 1.444 1.444
1.470 ID 1.440 1.444 1.444 1.444
1991 1995 1996 1997 2005
2.64 2.62 2.63 2.63 2.63
POL1l RPOLm TTMn AMOEBAo MAMEp
2.03 average: 0.990 3ID 2.03 average: 1.470 1.85 1.370 1.615 1.294 1.77 1.328 1.672 1.225 1.85 1.468 1.528 1.415
1990 1992 1999 2003 2005
2.56 2.82 2.58 2.48 2.64
TIP4P-FQq 1.85 0.820 TIP3P-FQ/Emr 1.85 0.889 TIP3P-FQ/Eqr 1.85 1.005
2.550 0.000 FQ 1.657 0.000 1.676 0.000
1994 2001 2001
2.29 2.33 2.39
POL5/TZs PPCt
1.494 1.060 FQ/ID 2001 1.01 0.00 FQ/COS 1996
2.56 2.50
1.85 1.320 2.14 0.66
a
The various models are characterized by their names, static dipole moments µ0x , diagonal elements Rll of the polarizability tensor, model class, year of publication, and average dipole moment 〈µx〉 as calculated from an average field 〈Ex(rO)〉 ) 1.62 V/Å at the oxygen atom (see Section 3). b From refs 7 and 8. c From ref 23. d From ref 24. e From ref 25. f From ref 22. g From ref 26. h From ref 27. i From ref 28. j From ref 29. k From ref 30. l From ref 31. m From ref 32. n From ref 33. o From ref 34. p From ref 35. q From ref 36. r From ref 37. s From ref 38. t From ref 39. u From ref 40.
the spring depends on the electric field E at the position rc of the opposite charge. The induced dipole is then given by µi ) qc(rO - rc) and varies linearly with the electric field: µi ) (qc2/ kc)E(rc). Thus, the polarizability R is isotropic and is given by qc2/kc. The models of the second class employ inducible point dipoles (IDs) µi ) RE with isotropic polarizabilities R. ID models including only one such dipole are essentially equivalent to the COS models (but less popular, because many MD codes can easily handle point charges while having problems with point dipoles). Three inducible dipoles (3IDs) are used in a set of computationally much more demanding approaches, in which the molecular polarizability becomes anisotropic. Here, only the model of Tsiper36 restricts the dipole-dipole interactions to different molecules, whereas all others32-35 also include intramolecular interactions. The third class contains models based on the fluctuating charge (FQ) approach. Within FQ, the values of the partial charges depend on the electrostatic potential Φ(ri) at the sites i of these charges. The FQ parameters are site specific electronegativities χ0i and “hardnesses” Ji. For details, we refer to the work of Rick et al.,37 but note that the location of the charges sites within the molecular x-y-plane restricts the FQ polarizability to this plane. The last two models listed in Table 1 are hybrids. The “polarizable point-charge” (PPC) model by Svishchev40 combines simplified versions of the COS and FQ schemes, whereas the so-called POL5/TZ model by Stern et al.39 combines ID and FQ concepts. Because of the common use of FQ, in these
Polarizability of Point-Polarizable Water Models hybrid models the polarizability is much larger within the molecular plane than in the perpendicular z-direction. Table 1 demonstrates that intense efforts have been spent on the development of polarizable models. Nevertheless, reviewing in 2002 the state of the art, Guillot came at that time to the following conclusion: “One has a taste of incompletion, if one considers that not a single water model available in the literature is able to reproduce with a great accuracy all the water properties. Despite many efforts to improve this situation very few significant progress can be asserted.” In his view, this apparent failure is due to “the extreme complexity of the water force field in the details, and their influence on the macroscopic properties of the condensed phase”. Thus, it is interesting to check whether the situation has meanwhile improved. For this check we will mainly concentrate on the size of the polarizability R, which is employed in the MM water models listed in Table 1, and on the size of the average dipole moment 〈µx〉, which is expected for these water models in a liquid water environment. 〈µx〉 depends on R, on the model-specific static dipole moment µ0x , and on the average field 〈E〉 acting on the water molecule in a liquid-phase environment. As will be shown by MM-MD simulations further below in Section 3, this field is oriented along the x-axis and has a value of approximately 〈Ex(rO)〉 ≈ 1.6 V/Å at the position rO of the oxygen atom. According to Table 1, the ID approaches unanimously employ for R and µ0x values close to the experimental ones and correspondingly arrive at values for 〈µx〉 near 2.63 D. For the complex 3ID approaches, one finds average dipole moments of 2.48 D and above. Here, the so-called POL1 model32 was characterized by a very small R, but a subsequent revision33 corrected this value to Rexp. Note that the recent “atomic multipole optimized energetics for biomolecular applications” (AMOEBA) potential35 gives the smallest value for 〈µx〉. Still further reduced (≈ 2.34 D) are the average dipole moments obtained for the FQ models. Interestingly, in the case of the COS models, the values for 〈µx〉 apparently decrease with time starting at 2.73 D in 1990 and recently approaching values below 2.4 D. This decrease is mainly due to a reduction of the value assumed for the polarizability R. This reduction came about by empirically fitting the model parameters to bulk-phase properties of liquid water and by abandoning the concept that R should be chosen close to Rexp. Among the recent COS models, for instance, the “simple water model (4-site) with negative Drude polarizability” (SWM4NDP) of Lamoureux et al.26 was optimized to reproduce the experimental values1,42-44 for the density, vaporization energy, self-diffusion constant, and static dielectric constant at room temperature and ambient pressure. According to the authors22 “the results clearly indicate that the value of R must be around 1.0 Å3 to yield reasonable liquid properties.” The resulting contradiction between an empirical choice for R, which is by 33% smaller than Rexp, demands an explanation. Here, Lamoureux et al. suggest22 that the reduced polarizability “is likely to arise form the energy cost of overlapping electronic clouds in the condensed phase due to Pauli’s exclusion principle opposing induction.” However, the polarizabilities calculated by quantum mechanics for the water dimer45 and water clusters46-48 are at most by 18% smaller than Rexp. Thus, the effects quoted by Lamoureux et al.22 cannot completely explain their empirical choice for R. In the context of an ongoing DFT/MM study on optimal water models of varying complexity, we have addressed the question of which value one should choose for the electronic polarizability R in polarizable MM water models. To obtain an answer,
J. Phys. Chem. B, Vol. 112, No. 19, 2008 6235 we had to distinguish changes of the molecular dipole moment, which result from field-induced variations of the average molecular geometry in the condensed phase, from changes that are caused by the interaction of the molecular electron system with a fluctuating external field. In this paper we concentrate on the latter effects, assuming that they are the dominant ones (which we have verified by test calculations). Therefore we will exclusively work with a fixed molecular geometry, for which we have chosen the experimental gas-phase finding. In the following section, we will present the computational methods applied by us in our search for a proper relation between the instantaneous polarization of a water molecule and an external electric field E(t). We first show how one can compute the electronic polarizability R of an isolated molecule, which is exposed to homogeneous electric fields. Furthermore, we show how one can derive R from DFT/MM calculations of a DFT molecule embedded in liquid water. Subsequently, we present our results and conclude the paper with a short summary and discussion. 2. Methods For our DFT and DFT/MM calculations of water molecules exposed to external electric fields, we employed the program packages CPMD49 and EGO/CPMD,50,51 respectively. Note that CPMD is a grid-based DFT code that uses plane wave basis sets for the expansion of the Kohn-Sham orbitals and pseudopotentials to model the core electrons. We applied the gradientcorrected exchange functional of Becke,52 the correlation functional of Perdew,53 and the norm-conserving pseudopotentials of Troullier and Martins54 to those water molecules, which represented the DFT fragments in our DFT/MM calculations. These DFT water molecules were fixed at the experimental gasphase geometry and were enclosed in a cubic box of edge length 9 Å, such that every nucleus was more than 3.7 Å distant from the faces of the box. The edges of the box were aligned with the coordinate axes shown in Figure 1. In CPMD, this box contains the grid for the plane wave basis set. The basis set was cutoff at 80 Ry. We denote this particular DFT approach as MT/BP. For comparison, a larger basis set as defined by a box dimension of 12 Å and a 120 Ry cutoff has also been considered. Because the computed dipole moments changed by only about 0.3%, we considered the 9 Å box with the 80 Ry cutoff as sufficiently accurate. The DFT/MM interface of the program package EGO/CPMD as described in Eichinger et al.50 provides a highly efficient way to import an external electrostatic potential into the DFT program CPMD. By this method, the potential ΦMM, which is generated by the MM fragment of a simulation system, is evaluated at the grid points rγ (γ ) 1, 2, ...) within the DFT box. Taking the gradient of the employed expressions,50 one immediately obtains the external electric field EMM(rγ) acting on the DFT fragment of the simulation system. The import of homogeneous electric fields into CPMD is particularly simple, (∞) because all those contributions ΦMM to the external potential, which are generated by MM charge distributions more distant than about 10 Å from a given DFT atom i, are represented by the coefficients of a second-order Taylor expansion50
1 (∞) t (∞) t (∞) Φ(∞) MM(rγ) ) Φi + (rγ - ri) Ei + (rγ - ri) Ti (rγ - ri) 2 Here, ri is the position of the DFT atom i closest to rγ, Φ(∞) is i the potential, E(∞) is the field, and T(∞) is the tensor of field i i gradients at ri. Thus, ignoring all DFT atoms but one, choosing
6236 J. Phys. Chem. B, Vol. 112, No. 19, 2008 zero for the corresponding coefficients Φ(∞) and T(∞) as well as a constant value for E(∞) generates a linearly varying external potential within the DFT box. To compute the components Rll, l ∈ {x, y, z} of the electronic polarizability tensor R for an isolated water molecule, we carried out a series of DFT calculations applying homogeneous fields of varying directions and strengths. Here, we chose electric fields E(∞) parallel to the axes of the Cartesian coordinate system shown in Figure 1. For the field components E(∞) l , we selected the range [-2, 2] V Å-1, because this range covers typical field strengths experienced by a water molecule in the bulk liquid at ambient temperature and pressure (cf. Section 3 below). For each field component l, the dipole moment µ was calculated (∞) by DFT at 19 different and equally spaced field values El,s (s ) 1, ..., 19) covering this range. Seventh order polynomial fits (∞) to the resulting dipole moments µ(El,s ) then yielded the 0 vacuum dipole moment µ and the respective components Rll of the polarizability tensor. Here, one of course has µ0 ) µ(0). Additionally, an isotropic polarizability R was calculated by a (∞) (∞) global fit to the absolute values |µi(|El,s |)| ) |µ(El,s ) - µ0| of the induced dipole moment. To answer the question, which value should be chosen for the polarizability R in MM-MD simulations of point-polarizable water models, we took a DFT/MM approach involving four steps: (i) MM-MD simulation systems of liquid water were set up and were equilibrated; (ii) a set of snapshots was extracted from these MM-MD trajectories covering an ensemble of environments representative for a solute water molecule embedded in liquid water; (iii) for each environment, the dipole moment of the solute molecule was calculated by DFT/MM; and (iv) the average electric field within the volume of that molecule was estimated. (i) Simulation Systems. Periodic orthorhombic dodecahedra with an inner radius of 21.8 Å were filled with 1970 TIP3P,19 TIP4P,19 and SPC/E20 water models, respectively. Using the EGO-MMII program package,51 these systems were equilibrated for 1 ns by MD in the NpT ensemble, i.e., at constant pressure p, temperature T, and number N of atoms. To keep the pressure at 1 × 10 5 Pa and the temperature at 293 K, we employed Berendsen thermo- and barostats55 at coupling times of 10 ps (T) and 0.1 ns (p), respectively. Additionally, the barostat was parametrized with an isothermal compressibility of 4.6 × 10-10 Pa-1, which is the experimental value for liquid water at the given conditions.1 In the MD simulations, the equations of motion were integrated by a multiple time-step algorithm56 with a basic time step of 1 fs. The long-range electrostatic interactions were treated by the structure-adapted fast multipole method57 combined with a reaction field correction, i.e., by the SAMM/ RF algorithm described in ref 51. The van der Waals interactions were cut off at 10 Å. The water molecules were kept rigid by the M-SHAKE algorithm58 with a relative tolerance of 10-6. (ii) Snapshots. From the last 100 ps of each equilibration trajectory, 20 snapshots were taken at a temporal distance of 5 ps. From each snapshot, 50 water molecules were randomly selected for the subsequent DFT/MM treatment, resulting in 1000 different environments for each of the three MM models summing up to a total of 3000 water molecules in different liquid surroundings s (s ) 1, ..., 3000). (iii) DFT/MM Dipole Moments. Choosing one of the thus selected water molecules as the DFT fragment and the remainder of the simulation system as the MM fragment, the dipole moment µs of the DFT fragment was calculated by the DFT/ MM approach outlined above. Here a slight technical difficulty arose in the treatment of the SPC/E model, because this model
Schropp and Tavan does not employ the experimental geometry. Therefore, the SPC/E molecules selected as DFT fragments had to be transformed into the experimental geometry first. For this purpose, we kept the position of the oxygen atom fixed and appropriately shifted the positions of the hydrogen atoms while retaining the molecular axis and plane. (iv) Average External Electric Fields. As explained further above, the EGO/CPMD interface provides simple access to the external electric field EMM(rγ) at the positions rγ of the grid points γ within the DFT box. Therefore, one can easily calculate average external electric fields 〈E〉 within the volume occupied by the electron density of the DFT fragment from these data. Here, however, the question arises how one should describe this volume. For the most simple modeling, we chose a Gaussian kernel of width σ and centered at the position rO of the oxygen atom. Denoting the number of grid points by G, for a given Gaussian width σ, the average external field then is G
〈E〉σ )
EMM(rγ) exp(-|rγ - rO|2/2σ2) ∑ γ)1 (1)
G
∑ exp(-|rγ |- rO| /2σ ) 2
2
γ)1
We evaluated 〈E〉σ for values of σ taken from the range [0.3 Å, 1.0 Å] in steps ∆σ ) 0.1 Å. Additionally, we evaluated 〈E〉σ for vanishing σ. Because limσf0〈E〉σ ) E(rO), this is simply the field at rO, which is directly accessible in MM-MD simulations. Once one has calculated for each of the 3000 snapshots s and each choice of σ the average fields 〈Es〉σ and, by DFT/ MM, the associated induced dipole moment µi,s ) µs - µ0 of the DFT fragment, one can plot these data for each σ into a correlation graph |µi(|〈E〉σ|)|. Linear regression subject to the constraint that µi ) 0 at 〈E〉σ ) 0 yields an isotropic polarizability Rσ through the slope of the straight line. Additionally, one obtains for each σ the standard deviation χσ, by which the induced dipole moments µi,s deviate from the regression line µi,σ ) Rσ〈E〉σ. 3. Results Employing the computational procedures described above, we first calculated by DFT/MM the dipole moment µ of a water molecule with fixed geometry in homogeneous electric fields El,s of varying directions and strengths. The zeroth and firstorder results of the polynomial fits to the resulting data µ(El,s) are the MT/BP predictions for vacuum dipole moment µ0 and for the diagonal components Rll of the electronic polarizability tensor (by symmetry, one has µ0 ≡ |µ0| ) µx0 and µy0 ) µz0 ) 0 for H2O). For the evaluation of the Rll, it is essential to use a fixed geometry in these DFT/MM calculations, because the Rll are experimentally determined through elastic and inelastic light scattering of water vapor,8 i.e., they represent purely electronic response properties. Table 2 compares the DFT results obtained for an isolated water molecule (exposed to homogeneous external fields) with corresponding experimental data. The table shows that our MT/ BP approach underestimates the observed vacuum dipole moment by 2.7% and overestimates the components of the polarizability tensor on the average by 7.5%. According to the MT/BP description, the isotropy of the polarizability R is even more pronounced than that reflected by the experimental data. But both sets of data indicate that the polarizability of a water molecule can be very well described by a single isotropic
Polarizability of Point-Polarizable Water Models
J. Phys. Chem. B, Vol. 112, No. 19, 2008 6237
TABLE 2: Comparison of DFT Results Obtained for an Isolated Water Molecule Exposed to Homogeneous Electric Fields with Experimental Dataa quantity
DFT (MT/BP)
experimentb
µ0 Rxx Ryy Rzz R
1.80 D 1.58 Å3 1.59 Å3 1.58 Å3 1.58 Å3
1.85 D 1.47 Å3 1.53 Å3 1.42 Å3 1.47 Å3
The results cover the vacuum dipole moment µ0, the diagonal elements Rll of the electronic polarizability tensor, and the isotropic polarizability R. See the text for further details and a discussion. b From refs 7 and 8. a
Figure 2. Cartesian components (crosses) of the dipole moment µ predicted by MT/BP for an isolated water molecule exposed to homogeneous fields directed parallel to the Cartesian coordinate axes l ∈ {x, y, z} (see Figure 1). The slope of the straight line indicates the isotropic polarizability Rhom. The gray bar marks the average electric field 〈E(rO)〉 at the oxygen atom as obtained by our MM-MD simulations for liquid water.
polarizability R. We call the MT/BP results for the isotropic polarizability Rhom and that for the vacuum dipole moment µ0hom. Figure 2 displays the dependence of the dipole moment µ on homogeneous fields, from which the data in Table 2 were derived, in greater detail and demonstrates that the assumptions of an isotropic polarizability and of linear response are valid within the considered range of field strengths. For all three field directions l, the components µl of the dipole moment are seen to vary linearly with the homogeneous external field for field strengths |El| smaller than about 1.5 V Å-1. For fields Ex > 0, which are directed parallel to the vacuum dipole moment µ0, the range of linear response extends up to field strengths of about 2 V Å-1, as demonstrated by Figure 2A. This case is important for the condensed phase, because here the water molecules will preferentially orient their permanent dipoles parallel to the external field. To indicate the typical field strengths to which a water molecule is exposed in the bulk liquid at room temperature and ambient pressure, we have evaluated the average electric field 〈E(rO)〉 at the oxygen atom from the 3000 MM-MD snapshots of liquid water described in Section 2. The values of the components 〈El(rO)〉 of this average field are indicated as gray bars in the three graphs. Whereas 〈Ey(rO)〉 and 〈Ez(rO)〉 vanish to a very good approximation, 〈Ex(rO)〉 is 1.62 V Å-1 underlining again that, in the liquid, the molecular dipoles are oriented parallel to the field. Furthermore, this value is well within the range of linear response. To estimate the electronic polarizability of water in the liquid phase, we calculated the dipole moments µs of an MT/BP water
Figure 3. Distribution of dipole moments µs of MT/BP water molecules calculated by DFT/MM for an ensemble of liquid MM water environments (s ) 1, ..., 3000) comprising equally sized sub-ensembles calculated with the TIP3P, TIP4P, and SPC/E potentials.
molecule exposed to 3000 different aqueous MM environments s as explained in Section 2. To ensure for this electronic property comparability with our gas-phase results, we have retained the fixed geometry, which we already employed in our DFT/MM calculations applying homogeneous external fields. Figure 3 shows a histogram of the resulting absolute values |µs| of the dipole moments. The average value of |µ| is 2.41 D, and the standard deviation is 0.11 D. Hence our DFT/MM result on the average dipole moment of a water molecule in a condensed phase environment is at the lower boundary of the range of values determined previously for this quantity (cf. our corresponding discussion in the Introduction). Furthermore, it is a little larger than the dipole moments of the surrounding MM water models, which range from 2.18 D (TIP4P) to 2.35 D (SPC/E, TIP3P). Using eq 1, we additionally evaluated for all environments the average electric fields 〈Es〉σ, which are generated by these environments, within Gaussian volumes of varying widths σ surrounding the oxygen atom of the DFT water molecule. Figure 4A shows the polarizability Rσ, which was derived by linear regression (cf. Section 2) from the correlation between the induced dipoles µi,s ≡ µs - µ0hom and the average external fields 〈Es〉σ, as a function of the widths σ of the Gaussian averaging volumes. According to the figure, Rσ increases with σ. The polarizability Rσ crosses the dotted line at σ ≈ 0.7 Å, which marks the value Rhom ) 1.58 Å obtained above from the response to homogeneous electric fields (cf. Table 2). At σ ) 0, i.e., if one evaluates the electric field E solely at the position rO of the oxygen atom within the DFT fragment, the polarizability assumes the value RO ) 1.08 Å3. This value is 32% smaller than the polarizability Rhom calculated by MT/BP for an isolated water molecule in homogeneous electric fields. The observed increase of the polarizability Rσ with the Gaussian width σ is caused by inhomogeneities of the electric field within the Gaussian averaging volumes. As demonstrated by Figure 4B, the ensemble average of the average electric fields 〈〈Ex,s〉σ〉s assumes large values at small or vanishing Gaussian volumes, where the field is solely averaged in a small vicinity of rO. Blowing up the Gaussian volume includes successively smaller field values in the averaging procedure such that 〈〈Ex,s〉σ〉s decreases with σ. As we will show below, this trend also holds for the individual members 〈Es〉σ ≈ µi,s/Rσ of the MM snapshot ensemble. The approximate inverse proportionality between 〈Es〉σ and Rσ then proves the above claim. Figure 4C illustrates at which quality the induced DFT/MM dipoles µi,s can be represented by the linear response expression Rσ〈Es〉σ. The figure shows the σ-dependence of the standard deviation χσ between the induced dipoles µi,s and the values predicted by Rσ〈Es〉σ. Large values of χσ indicate poor approximations, and smaller values point to improved descriptions. According to the graph, χσ becomes minimal, and the linear response relation becomes best at σopt ) 0.7 Å, where χσ
6238 J. Phys. Chem. B, Vol. 112, No. 19, 2008
Schropp and Tavan
Figure 6. Correlation graph |µi(|E(rO)|)| for the DFT/MM dipole moments and the corresponding fields at rO. The gray line depicts a linear regression, which was constrained to pass |µi| ) 0 at |E(rO)| ) 0 and was employed to determine RO ) 1.08 Å3 from the slope. The dashed gray line represents an unconstrained linear regression to the data, which yields µˆ 0i ) 0.18 D and Rˆ O ) 0.78 Å3 as parameters. For a discussion, see the text.
Figure 4. (A) The DFT/MM polarizability Rσ as a function of the Gaussian width σ. The polarizabilities Rσ were obtained by linear regressions from the data sets {(µi,s, 〈Es〉σ)|s ) 1, ..., 3000} correlating the induced dipole moments with electric fields averaged over Gaussian molecular volumes of widths σ. The dotted line indicates the polarizability Rhom ) 1.58 Å calculated by MT/BP for varying homogeneous fields. (B) σ-dependence of the ensemble average 〈...〉s of the average electric fields 〈Ex,s〉σ within the Gaussian volumes of widths σ. (C) σ-dependence of the standard deviations χσ between the induced dipoles µi,s calculated by DFT/MM and the regression lines µi,σ ) Rσ〈E〉σ.
Figure 5. Correlation graph |µi (|〈E〉σopt|)|. The gray line depicts the linear regression employed to determine Rσopt from the slope. The dashed light-gray line represents the linear response |µi| ) Rhom|〈E〉σopt|.
assumes a value of 0.02 D (note that this value for σopt represents only a rough estimate, because we have sampled the σ values in quite coarse steps of ∆σ ) 0.1 Å). Interestingly, the thus determined value σopt is very close to the point at which Rσ is just equal to Rhom (cf. Figure 4A). Thus the DFT/MM results on the dipole moments of MT/BP water molecules in MM water can be nicely reproduced by the simple linear response relation Rhom 〈Es〉σopt, in which one uses the value Rhom for the polarizability, which was derived for an isolated MT/BP water molecule exposed to homogeneous electric fields. Note that this quantitative comparison of polarizabilities has been enabled by the use of identical and fixed geometries in the condensed phase and homogeneous external field calculations, respectively. For the specific choice σ ) σopt, the correlation graph |µi(|〈E〉σopt|)| in Figure 5 demonstrates the quality at which the induced DFT/MM dipole moments µi,s represent a linear function of the average electric fields 〈Es〉σopt. The optimal linear
relation is given by the slope Rσopt ) 1.53 Å3 of the linear regression (gray line), but the linear response relation based on Rhom ) 1.58 Å3 (dashed light-gray line) approximates the data also very well. For the latter linear response relation Rhom 〈Es〉σopt, the relative error in the description of the DFT/MM ensemble of induced dipole moments µi,s is about 5% (for Rσopt it is 3%). For other values of σ, the relative error becomes much larger when using the linear response relation Rhom 〈Es〉σ to estimate µi,s. At σ ) 0, for instance, i.e., when measuring the electric field exclusively at rO, the relative error is 50%. This huge difference of relative errors immediately explains why one necessarily gets suboptimal descriptions of the bulkphase properties in MM-MD simulations with point-polarizable water models, if one employs the experimental polarizability for the computation of µi through the linear response Rexp E(rO) to the field at the oxygen atom. For an optimal description, one should instead employ the linear response Rexp 〈Es〉σopt to the average field around rO, which, however, is essentially inaccessible in such MM-MD simulations. On the other hand, Figure 4B has shown that the ensemble average |〈〈Es〉σopt〉s| is by a certain factor (γ ) 0.724) smaller than |〈E(rO)〉s|, such that one should get a close to optimal description by using the linear response
µi ) γ Rexp E(rO)
(2)
Because of our coarse sampling of the σ axis, the value σopt and, consequently, the scaling factor γ are not very precisely determined through the quotient of ensemble averaged fields. A more precise estimate for γ can be obtained through the quotient RO/Rhom, where the polarizability RO is calculated from the constrained linear regression at σ ) 0. The reasoning behind this approach is that, in our DFT description, Rexp is replaced by Rhom and γRhom corresponds to RO. The resulting representation of γ then is independent of the degree of effort spent to accurately determine σopt. Figure 6 shows the correlation graph |µi(|E(rO)|)| for electric fields evaluated at rO. The constrained linear regression (gray line) yields the value RO ) 1.08 Å3 quoted further above in the discussion of Figure 4A, which is required for the computation of the scaling factor γ. Thus, we find for the scaling factor γ the value
γ≡
RO ) 0.684 Rhom
(3)
Polarizability of Point-Polarizable Water Models
Figure 7. Correlation between the average electric fields |〈Es〉σopt| and the electric fields |EO,s| ≡ |Es(rO)| measured at rO in the ensemble of MM environments s. The gray line represents a linear regression. For a discussion, see the text.
In Figure 6 the induced dipole moments |µi,s| are seen to scatter much more strongly around the constrained linear regression than in Figure 5, illustrating the difference between the standard deviations χ0 ) 0.07 D and χσopt ) 0.02 D (cf. Figure 4C). The reason for this increased scatter becomes apparent by inspection of Figure 7, which correlates the volume average fields |〈Es〉σopt| with the spot checks |EO,s| at rO. Quite obviously, the spot checks |EO,s| are systematically larger than the volume averages |〈Es〉σopt|. According to the results that are displayed in Figure 5 above, the volume averages closely represent the fields sensed by a quantum mechanically described water molecule. Thus, the spot checks |EO,s| do not only overestimate these fields by neglecting the inhomogeneity of the field within the volume occupied by a water molecule, but additionally neglect the fluctuations of the field inhomogeneity around the average inhomogeneity. In Figure 7, the neglect of the average field inhomogeneity is documented by the slope of the regression line, which is 1.38 instead of one, and the neglect of the fluctuations shows up through the scatter around that line. Whereas the systematic overestimate can be corrected using the scaling eq 2, the errors caused by the fluctuations remain. However, despite the neglect of the fluctuations of the field inhomogeneity, the linear response RO Es(rO) still describes the induced DFT/MM dipoles µi,s quite well. The scatter of the data in Figure 6 can be quantified as a relative error of 11%, which is a little larger than the 3% error found in connection with Figure 5 for the linear response to the average field at σopt. As indicated by the unconstrained regression line (dashed gray line) in Figure 6, one may possibly obtain an improved description of the electronic polarizability in point-polarizable MM water models by adding a constant off-set
∆µ0 ≡ µˆ 0i
Rexp ) 0.17 D Rhom
to the experimental vacuum dipole moment µ0exp and by using an even smaller scaling factor
γˆ ≡
Rˆ O ) 0.493 Rhom
to reduce the experimental polarizability. However, because this modified approach does not yield the experimental vacuum dipole moment for the gas phase, it seems to be restricted to the specific case of bulk liquid water and apparently lacks transferability to other environments. Therefore we suggest that one should employ the combination of µ0exp with the scaled linear response relation defined by eqs 2 and 3 in pointpolarizable MM models for water, which are designed to be transferable.
J. Phys. Chem. B, Vol. 112, No. 19, 2008 6239 The value determined for the scaling factor γ, by which the experimental polarizability Rexp should be corrected in pointpolarizable MM water models, depends, of course, on the MM water models employed in the DFT/MM descriptions. In our study we have chosen three different MM models (TIP3P, TIP4P, and SPC/E) to reduce the model bias. The size of that bias can be estimated by evaluating γ through eq 3 individually for the sub-ensembles associated with these three models. For TIP3P, we find a 4.5% larger γ. For TIP4P and SPC/E, the values are 0.9% and 1.9% smaller than our average. Thus, the modifications to γ specific for a given MM model are small. Furthermore, the value of γ solely accounts of the field inhomogeneity within the volume of a given molecule and does not cover further effects, which additionally can modify the polarizability of water molecules in the condensed phase. For instance, the average geometry of a water molecule will change upon transfer from the gas phase into the condensed phase, and the electronic polarizability will be slightly affected by this change of geometry. Similarly, the altered average geometry will be associated with a changed static dipole moment. Preliminary estimates from DFT/MM dynamics simulations of a single DFT water molecule in an MM water environment point toward a geometry-mediated increase of the polarizability by about 3% and of the static dipole moment by about 1.2%. Finally, many-body effects such as charge transfer or Pauli repulsion within water clusters may cause further modifications. Having identified the field inhomogeneity as the key effect, these issues can now be systematically addressed in subsequent studies. 4. Summary and Discussion Studying by DFT/MM methods the electrostatic response of a stiff H2O molecule to external electric fields, we found that this response is well described by the polarizability Rhom, which we calculated for an isolated water molecule exposed to homogeneous electric fields. For the liquid phase, in which every water molecule is exposed to strongly inhomogeneous and varying external fields, we found that the electronic polarizations of the water molecules represent responses to fields 〈E〉, which are aVerages over the volumes of these molecules. Because of the field inhomogeneities, these average fields 〈E〉 are sizably smaller than the fields E(rO) at the oxygen atoms, which are readily accessible in MM-MD simulations (as opposed to the average fields 〈E〉) and, therefore, were previously used in MMMD simulations with point-polarizable water models. Because computational efficiency dictates that one has to employ E(rO) in such MM-MD simulations, one cannot use the experimental polarizability Rexp, but must use a reduced “effective” polarizability Reff ) γ Rexp. Here, the required scaling factor γ is obtained by eq 3 from Rhom and from the correlation between the induced dipoles µi,s and fields Es(rO) calculated by DFT/MM for an ensemble of liquid-phase environments s. With our specific choice of the MT/BP method for the DFT calculations and of the TIP3P, TIP4P, and SPC/E potentials for the liquid-phase MM environments, we found for γ a value of 0.684. On the basis of our results, we thus suggest that one should combine the static dipole moment µ0exp ) 1.855 D with the effective polarizability Reff ) 1.005 Å3 in the construction of transferable point-polarizable water models. As also outlined in Section 3, one could also consider for water models specific to liquid water the combination µ0 ) 2.024 D with Reff ) 0.725 Å3. Remarkably, our DFT/MM-based suggestion for the electrostatic properties (µ0, Reff) of a transferable water model closely
6240 J. Phys. Chem. B, Vol. 112, No. 19, 2008 agrees with the completely empirical optimization results of Lamoureux et al.22,26 These authors also chose the experimental dipole moment µ0exp and combined it with polarizabilities Reff deviating by only 3.7% (2003) and 2.7% (2006) from our result (cf. Table 1). Our analysis has clearly demonstrated that the scaling factor γ, which reduces Rexp to Reff, exclusively accounts for the field inhomogeneities in the liquid phase, which were not considered in previous analyses. These inhomogeneities have the particular property that the external field acting on a water molecule in the liquid assumes its maximum at the oxygen atom. It is as if the electrostatics of the water molecules would steer the oxygen atoms to locations of maximal field strengths. Therefore, measuring the field strengths at these locations is a particularly bad choice for the evaluation of the induced dipole moments by linear response with an unscaled polarizability, e.g., with Rexp. As a final result, we thus can entirely explain the strongly scaled Reff suggested by Lamoureux et al.22 through the neglected field inhomogeneities without having to refer to effects originating from “Pauli’s exclusion principle opposing induction”, which anyway cannot account for the whole reduction of R (cf. Section 1). Acknowledgment. This work was supported by the Deutsche Forschungsgemeinschaft (SFB 533/C3). References and Notes (1) Eisenberg, D.; Kauzmann, W. The Structure and Properties of Water; Clarendon Press: Oxford, 1969. (2) Dougherty, R. C. J. Chem. Phys. 1998, 109, 7379-7393. (3) Errington, J. R.; Debenedetti, P. G. Nature 2001, 409, 318-321. (4) Guillot, B. J. Mol. Liq. 2002, 101, 219-260. (5) Head-Gordon, T.; Hura, G. Chem. ReV. 2002, 102, 2651-2669. (6) Benedict, W. S.; Gailar, N.; Plyler, E. K. J. Chem. Phys. 1956, 24, 1139-1165. (7) Clough, S. A.; Beers, Y.; Klein, G. P.; Rothman, L. S. J. Chem. Phys. 1973, 59, 2254-2259. (8) Murphy, W. F. J. Chem. Phys. 1977, 67, 5877-5882. (9) Gregory, J. K.; Clary, D. C.; Liu, K.; Brown, M. G.; Saykally, R. J. Science 1997, 275, 814-817. (10) Sprik, M. J. Chem. Phys. 1991, 95, 6762-6769. (11) Soetens, J.-C.; Costa, M. T. C. M.; Millot, C. Mol. Phys. 1998, 94, 577-579. (12) Silvestrelli, P. L.; Bernasconi, M.; Parrinello, M. Chem. Phys. Lett. 1997, 277, 478-482. (13) Silvestrelli, P. L.; Parrinello, M. J. Chem. Phys. 1999, 111, 35723580. (14) Badyal, Y. S.; Saboungi, M.-L.; Price, D. L.; Shastri, S. D.; Haeffner, D. R.; Soper, A. K. J. Chem. Phys. 2000, 112, 14559-14560. (15) Gubskaya, A. V.; Kusalik, P. G. J. Chem. Phys. 2002, 117, 52905302. (16) Hohenberg, P.; Kohn, W. Phys. ReV. 1964, 136, B864-B870. (17) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133-A1138. (18) Tavan, P.; Carstens, H.; Mathias, G. Molecular dynamics simulations of proteins and peptides: Problems, achievements. and perspectives. In Protein Folding Handbook; Buchner, J., Kiefhaber, T., Eds.; WileyVCH: Weinheim, Germany, 2005; Vol. 1.
Schropp and Tavan (19) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (20) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269-6271. (21) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910-8922. (22) Lamoureux, G.; MacKerell, A.; Roux, B. J. Chem. Phys. 2003, 119, 5185-5197. (23) Straatsma, T. P.; McCammon, J. A. Chem. Phys. Lett. 1990, 167, 252-254. (24) Yu, H.; Hansson, T.; van Gunsteren, W. F. J. Chem. Phys. 2002, 118, 221-234. (25) Yu, H.; van Gunsteren, W. F. J. Chem. Phys. 2004, 121, 95499564. (26) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D. Chem. Phys. Lett. 2006, 418, 245-249. (27) Kozack, R. E.; Jordan, P. C. J. Chem. Phys. 1991, 96, 3120-3130. (28) Brodholt, J.; Sampoli, M.; Vallauri, R. Mol. Phys. 1995, 86, 149158. (29) Chialvo, A. A.; Cummings, P. T. J. Chem. Phys. 1996, 105, 82748281. (30) Dang, L. X.; Chang, T.-M. J. Chem. Phys. 1997, 106, 8149-8159. (31) Paricaud, P.; Predota, M.; Chivalvo, A. A.; Cummings, P. T. J. Chem. Phys. 2005, 122, 244511-244514. (32) Caldwell, J.; Dang, L. X.; Kollman, P. A. J. Am. Chem. Soc. 1990, 112, 9144-9147. (33) Dang, L. X. J. Chem. Phys. 1992, 97, 2659-2660. (34) Burnham, C. J.; Li, J.; Xantheas, S. S.; Leslie, M. J. Chem. Phys. 1999, 110, 4566-4581. (35) Ren, P.; Ponder, J. W. J. Phys. Chem. 2003, 107, 5933-5947. (36) Tsiper, E. V. Phys. ReV. Lett. 2005, 94, 013204-4. (37) Rick, S.; Stuart, S.; Berne, B. J. Chem. Phys. 1994, 101, 61416156. (38) Ando, K. J. Chem. Phys. 2001, 115, 5228-5237. (39) Stern, H.; Rittner, F.; Berne, B.; Friesner, R. J. Chem. Phys. 2001, 115, 2237-2251. (40) Svishchev, I. M.; Kusalik, P. G.; Wang, J.; Boyd, R. J. J. Chem. Phys. 1996, 105, 4742-4750. (41) Yu, H.; van Gunsteren, W. F. Comput. Phys. Commun. 2005, 172, 69-85. (42) Jancso, G.; Van Hook, W. Chem. ReV. 1974, 74, 689-750. (43) Krynicki, K.; Green, C.; Sawyer, D. Faraday Discuss. Chem. Soc. 1978, 66, 199-208. (44) Ferna´ndez, D.; Mulev, Y.; Goodwin, A.; Sengers, J. J. Phys. Chem. Ref. Data 1995, 24, 33-70. (45) in het Panhuis, M.; Popelier, P. L. A.; Munn, R. W. J. Chem. Phys. 2001, 114, 7951-7961. (46) Morita, A.; Kato, S. J. Chem. Phys. 1999, 110, 11987-11998. (47) Morita, A. J. Comput. Chem. 2002, 23, 1466-1471. (48) Tu, Y.; Laaksonen, A. Chem. Phys. Lett. 2000, 329, 283-288. (49) CPMD V3.9; Copyright IBM Corp. 1990-2004, Copyright MPI fu¨r Festko¨rperforschung Stuttgart 1997-2001; see also www.cpmd.org. (50) Eichinger, M.; Tavan, P.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1999, 110, 10452-10467. (51) Mathias, G.; Egwolf, B.; Nonella, M.; Tavan, P. J. Chem. Phys. 2003, 118, 10847-10860. (52) Becke, A. D. Phys. ReV. A 1988, 38, 3098-3100. (53) Perdew, J.; Yue, W. Phys. ReV. B 1986, 33, 8800-8802. (54) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993-2005. (55) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (56) Eichinger, M.; Grubmu¨ller, H.; Heller, H.; Tavan, P. J. Comput. Chem. 1997, 18, 1729-1749. (57) Niedermeier, C.; Tavan, P. J. Chem. Phys. 1994, 101, 734-748. (58) Kra¨utler, V.; van Gunsteren, W. F.; Hu¨nenberger, P. J. Comput. Chem. 2001, 22, 501-508.
