Transferable ab Initio Dipole Moment for Water: Three Applications to

Oct 5, 2015 - The average monomer dipole moments for 300 K water and 0 K ice Ih are 2.94 ... at 300 K from 0 to 4000 cm Is in Near-Quantitative Agreem...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCB

Transferable ab Initio Dipole Moment for Water: Three Applications to Bulk Water Hanchao Liu, Yimin Wang, and Joel M. Bowman* Cherry L. Emerson Center for Scientific Computation and Department of Chemistry, Emory University, Atlanta, Georgia 30322, United States ABSTRACT: We recently reported a second-generation, ab initio dipole moment surface (DMS) for water and applied it successfully to the IR spectrum of liquid water at 300 K. Here the transferability of this DMS is demonstrated in three applications. One is the distribution of monomer dipole moments, considering two definitions, and effective atomic charges of liquid water at 300 K and also for a model of ice Ih at 0 K. The second one is a calculation of the dielectric constant of liquid water at 280, 300, 320, 340, 360 K, and the third one is correcting the intensities of the IR spectrum of liquid water at 300 K, obtained using the q-TIP4P/F potential, bringing them into much improved agreement with experiment. For the purpose of obtaining statistical ensembles we use molecular dynamics simulations with the TIP4P+E3B water model developed by Skinner and co-workers. The average monomer dipole moments for 300 K water and 0 K ice Ih are 2.94 and 3.54 D, respectively, in good agreement with literature values. Effective monomer charge distributions are derived from the monomer dipoles and give average values of −1.02 e for O and 0.51 e for H in liquid water, which are also in agreement with values reported from experiment. The calculated dielectric constant of liquid water at the above five temperatures is compared to experiment and is roughly 10−15% lower than experiment.

1. INTRODUCTION Developing accurate models for the potential and electrostatic properties of water is an active and widespread research field. An important example of a property is the dipole moment, since it is directly related to the intensity of the IR spectrum and other properties such as the dielectric constant. Most models for describing the dipole moment are based on fixed charges, e.g., rigid monomer TIP4P,1 SPC2 and flexible versions, q-TIP4P/F,3 and q-SPC/wF.4 More sophisticated treatments of the electrostatics include induction effects by accounting for the polarizability of the water monomers. There are many variations of polarizable models, e.g., the TTMn family of potentials,5 and recent versions of the AMOEBA potential.6 It is worth stressing that, in general, the dipole and induced dipoles are an integral component of the potential in these models. However, the efficient fixed-charge dipole models do not capture the significant change in the IR intensities of water clusters and liquid water compared to the isolated gasphase spectrum in the bend and OH-stretch bands. This is wellknown and was noted previously.3 In fixed-charge water models, charges are placed on several sites of a water molecule. Usually, the number of sites is 3 or 4 but can be extended to as many as 6. The charges in a model are parameters that are determined to fit selected experimental properties of bulk water/ice. In the three-site water model such as SPC, a negative charge, usually of the order of −0.8 e, is on the oxygen and a positive charge is on each hydrogen to neutralize the molecule. In the four-site water model such as TIP4P, while a positive charge is on each hydrogen, a negative charge is placed on a virtual site M, which is on the bisector of © XXXX American Chemical Society

the bond angle. The magnitude of the negative charge is usually between −1 and −1.2 e. Due to the parametrization to bulk properties, the monomer dipole moment from these models is usually in the range of 2.2−2.6 D, significantly larger than the gas phase dipole moment of 1.85 D. In these models, the total dipole moment for a finite water cluster is given by the vector sum of the monomer dipoles. One “first-principles” way to describe the dipole moment of condensed-phase water is to use a many-body representation. In this approach, the exact dipole moment is represented by a sum of one-, two-body, etc. dipole moments. In this representation, analogous to the one for the potential,7 the one-body dipole moment is the one for the isolated monomer, and the two-body one, termed the intrinsic two-body dipole moment, is the difference of the total dipole moment of the water dimer minus the two monomer dipole moments, etc. The flexible dipole moment surface we developed is truncated at one-body and two-body level. To the best of our knowledge, this representation of the flexible ab initio dipole moment was first proposed by Wang et al.7 It has been shown in tests to reproduce the dipole moment of water clusters up to the 20mer to within 10% or less. It is perhaps worth noting that the one-body dipole, unlike the empirical monomer dipole in fixedcharge models, is the dipole moment of an isolated water Special Issue: Bruce C. Garrett Festschrift Received: September 21, 2015 Revised: October 1, 2015

A

DOI: 10.1021/acs.jpcb.5b09213 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

moments). The recent study by Skinner and co-workers17 reported that ϵ is 52 from TIP4P model and is 56 from their rigid three-body empirical potential, E3B, addition to TIP4P. Clearly, the results are significantly below the experimental value of 78.5, at 298 K. They pointed out these models do not include polarization effects which cause the underestimation of the ϵ. They proposed a linear model to describe the induced dipole as a function of electric field, and the coefficient is empirically determined to reproduce the experimental ϵ. ϵ decreases as temperature increases. The dielectric constant for ice Ih was also reported in this paper; it is greater than the values of liquid water and behaves similarly with temperature. The third application of the dipole moment surface is to the IR spectrum of liquid water, with a focus on the relative intensity of the OH-stretch and monomer bend bands. As noted above, we used the DMS, along with the ab initio WHBB PES and the local-monomer model, to calculate the IR spectrum of liquid water at 300 K.13 The agreement with experiment is excellent for both the frequency and intensity. Here we use the DMS but with the q-TIP4P/F potential to calculate the IR spectrum, again using the local-monomer model. The q-TIP4P/F potential is an efficient, flexible version of TIP4P potential, which belongs to the fixed-charge class of water models.3 It is computationally efficient and is able to predict reasonably accurate vibrational frequencies of liquid water; however, it is unable to give the correct intensity.3,13 We are motivated to improve the intensity by hybridizing the qTIP4P/F potential with the ab initio dipole moment surfaces. In these three applications (with the exception of the one to ice Ih) NVT molecular dynamics (MD) calculations are done to obtain the water monomer ensembles at a given temperature. The TIP4P+E3B potential is used in these MD calculations. Then these ensembles are used in subsequent calculations in the three applications, using the ab initio DMS. This is explained in detail below. The paper is organized as follows. In section 2, details of the MD calculations are given, followed by details of the ab initio dipole moment surface, the monomer dipole moment, dielectric constant, and local-monomer IR calculations. In section 3, the results and discussion are given for each application. A summary and concluding remarks are given in section 4.

molecule. In the latest version, generation 2 of our DMS, we use a spectroscopically accurate dipole moment surface for the water molecule, developed by Tennyson and co-workers.8 The intrinsic two-body DMS is a fit to roughly 30 000 ab initio data points using a fitting basis of polynomials (in Morse variables) which ensures that the dipole moment transforms covariantly under permutation of same-type atoms. Because the intrinsic two-body dipole is directly derived from ab initio calculations, it includes the effect of polarization and intermolecular charge transfer. Previously, the DMS has been successfully and almost exclusively applied to the infrared spectra of water clusters on small water clusters,9−11 ice,11,12 and liquid water,13 along with the ab initio WHBB water potential energy surface7 and the quantum local-monomer model.9 We note that recently another many-body flexible DMS using a similar approach has been reported by Medders and Paesani, which also includes higher-order polarization contributions to the dipole moment.14 Clearly, the WHBB potential (which is also represented by a many-body expansion) and the one-, two-body dipole moment are completely uncoupled from each other. Thus, they are independently transferable; e.g., the dipole moment surface can be used with any potential for water, for example, a fixed-charge one to obtain the IR spectrum. The dipole moment surface can also be used to obtain electrical properties of water clusters of arbitrary size including the bulk. For example, configurations of bulk water simulations can also be used to obtain average electrical quantities such as the total dipole moment, model monomer dipoles, the dielectric constant, etc. In this paper, we present three applications using the dipole moment surface. First, we consider the distribution of model monomer dipole moments in liquid water and ice Ih. This property is very useful because it reflects the change of electrical environment from gas phase to condensed phase. The definition of monomer dipole moment is subject to the models and methods used to calculate it. Although the monomer dipole lacks a universal definition, there is consensus that the monomer dipole moment increases from the value in gas phase, 1.85 D, to around 2.5−3.0 D.6,7,15−19 The increase is an indication of solvent effects, such as polarization and charge transfer. X-ray diffraction experiments provide an indirect determination of the monomer dipole of 2.9 ± 0.6 D by fitting a charge transfer parameter of 0.5 e along each OH bond.15 A number of theoretical studies obtain a value for the dipole in the range 2.4−3.1 D.6,16,17 Wang et al. proposed a simple approach to obtain an average monomer dipole moment from the total dipole moment (which is a rigorous quantity) and applied it to clusters.7 Kemp and Gordon determined the enhancement of the monomer dipole in bulk water from the gas phase by using localized charge distribution methods.18 An average monomer dipole moment of 3.0 D was reported by Silvestrelli and Parrinello, using a localized orbital analysis in 318 K ab initio molecular dynamics.19 Recently, Shi et al.17 reported the distribution of monomer dipoles for liquid water at 303 K and ice Ih at 245 K and also the dielectric constant of ice and liquid water with the introduction of a scaling parameter to match experiment. This work will be discussed in more detail below. The second application is a calculation of the dielectric constant, ϵ, of liquid water. Water has a large dielectric constant, about 78 at 298 K. This large dielectric constant is important for stabilizing ions in solutions. This property is also one of the standards for assessing water potentials (and dipole

2. METHODS AND CALCULATION DETAILS 2.1. MD Simulation Details of Liquid Water and the Model of Ice Ih. Canonical ensembles of water monomers for liquid water are obtained by classical NVT trajectories of 500 water molecules. The water potential for these calculations is the TIP4P two-body potential plus the semiempirical rigid three-body E3B potential.20,21 The simulation is performed at five temperatures, 280, 300, 320, 340, 360 K. The Nose− Hoover thermostat is used to control the temperature, with a time constant of 5 ps.22,23 The usual periodic boundary conditions are applied. The box size is adjusted to reproduce the experimental density at each temperature. After a 2 ns equilibration phase, a 10 ns production trajectory is collected at each temperature. The interval between two adjacent frames is 100 fs. A modified version of GROMACS 4.5.5, with the E3B water model included, provided by the Skinner’s group, is used for all the MD simulations.20,21 The model for ice Ih at 0 K is adopted from our previous applications.12 This model contains 192 water molecules, with a structure proposed by Hayward and Reimers.24 The “proton” B

DOI: 10.1021/acs.jpcb.5b09213 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (that is, the H atom) is randomly disordered with a constraint to ensure zero total dipole moment. The original structure is further optimized on the WHBB water potential. We have identified 105 monomers in the core region of this model and 87 monomers on the surfaces.12 2.2. Ab Initio Dipole Moment Surface. A detailed description of our latest dipole moment surface is described elsewhere.13 In brief, the DMS is a many-body representation using one-body and two-body dipoles. Thus, the total dipole moment for N monomers, M⃗ , is given by N

M⃗ =

Clearly, this definition results in a smaller monomer dipole moment. Results from both definitions will be presented and discussed in section 3. The calculation of monomer dipoles is done at 300 K. The distribution of monomer dipole moment magnitudes is obtained straightforwardly from the MD ensemble of dipole moments of each monomer. Specifically, we select 100 MD snapshots separated by 100 ps and then select 100 core monomers in each snapshot. Therefore, a total of 10 000 monomer configurations are sampled to obtain the monomer dipole distribution. A Gaussian-smoothed distribution is obtained by replacing each data point with a Gaussian function, centered at the magnitude of each monomer dipole moment and with a fwhm of 0.1 D. The final distribution is obtained by adding all the Gaussian-broadened data. Note that choosing a fwhm of 0.05 D produces a distribution that is visually nearly identical to the one using 0.1 D. Note that in the TIP4P+E3B MD simulation, water molecules are rigid, using the gas-phase equilibrium geometry of a bond length of 0.957 Å and a bond angle of 104.5°. A more realistic geometry, as pointed out previously,25 is the ZPEaveraged one. Thus, the OH bond length is 0.9716 Å, and the bond angle is 104.7°.25 This is implemented trivially using this configuration for every monomer in each snapshot with this configuration, with the monomer’s center of mass kept unchanged. This geometry is used in the calculation of monomer and total dipole moments. Finally, it is of interest to determine atomic charges on each monomer that reproduces the monomer dipole moment. This can be done easily for H2O once the dipole moment is determined. Specifically, the charges of the three atoms are calculated by solving the linear equation

N

∑ μ1⃗ ‐body (i) +

∑ μ2⃗ ‐body (i ,j)

i=1

i>j

(1)

For the one-body dipole, we use the spectroscopically accurate DMS developed by Tennyson and co-workers, denoted as LTP2011.8 The two-body dipole is a permutationally covariant fit of roughly 30 000 MP2/aug-cc-pVTZ dipole moments. The fitting basis is a set of polynomials in Morse variables up to fifth order. Here, we report another fit of the two-body dipole using the same data set but with a smaller fitting basis, i.e., a set of polynomials up to fourth order. The purpose of this new fourth-order DMS is to increase the computational efficiency. It will be applied, along with the more precise fifth-order fit, to the calculation of the IR spectrum of liquid water using the qTIP4P/F potential. Details of these calculations are given in section 2.5. 2.3. Monomer Dipole Moment and Monomer Effective Charges. As mentioned in the Introduction, the monomer dipole moment, while conceptually useful, does not have a rigorous definition. Within the many-body representation of the total dipole, two physically reasonable definitions of the monomer dipole are proposed. In the first definition, the dipole moment of monomer m is written as → ⎯μ = μ (m ) + m 1⃗ ‐body

∑ μ2⃗ ‐body (m ,k) k≠m

X mqm⃗ = μm⃗

in which the Xm is a 3 × 3 matrix of the monomer’s Cartesian coordinates. qm⃗ is a 3 × 1 vector of the three atomic charges of the monomer. μ⃗m is a 3 × 1 vector of the three components of the monomer dipole. Note that this equation does not constrain the net charge of a monomer to be zero, and we comment on this below. Also, it is worth noting that these atomic charges are not the same as the “charges” that are determined by the covariant least-squares fitting of the dipole moment.7,13 The fitting is done to the dipole moment, which in the case of the two-body dipole has six “charges” and of course three dipole components. The fitted uses the constraint that the sum of the fitting “charges” is zero. This is not sufficient to uniquely determine the six “charges”. They generally come out reasonable, and in the limit of separated monomers they do agree with the charges determined independently for the monomer.8 In addition, the fitted charges of the intrinsic twobody dipole go to zero as the dimer separates.13 2.4. Dielectric Constant. The dielectric constant is calculated using the usual expression,

(2)

in which the first term is the one-body dipole for the monomer and the second term is the sum of the two-body terms involving monomer m. This definition is physically intuitive, as the one-body dipole reflects the electron density distribution of a water molecule by itself and all the intrinsic two-body dipoles reflect the change of monomer’s electron density distribution caused by the presence of other monomers. As mentioned earlier, the intrinsic two-body dipole is obtained directly from fitting ab initio calculations, so it implicitly describes effects such as polarization and intermolecular charge transfer. Note that with the above definition of the monomer dipole moment, the sum of all monomer dipole does not equal the total dipole of the cluster given by eq 1. This is because summing the monomer dipoles, using the above definition, double-counts the two-body terms. (This is not an issue in fixed-charge models because there is only effective one-body dipole and the sum of those does give the total dipole.) If we wish to satisfy the condition that the sum of monomer dipole equals the total dipole moment, then this is achieved using the modified, second definition below. 1 → ⎯μ = μ (m ) + m 1⃗ ‐body 2

∑ μ2⃗ ‐body (m ,k) k≠m

(4)

ϵ=1+

4π (⟨M2⟩ − ⟨M⃗ ⟩2 ) 3Vk bT

(5)

in which kb is the Boltzmann constant, T is the temperature, and V is the volume of the box and is temperature dependent. The quantity in the parentheses is the variance of the total dipole moment. M⃗ is the total dipole from the DMS, calculated from eq 1. M2 is defined as M2 = M⃗ ·M⃗ . Note that no periodic boundary conditions are applied when computing the total

(3) C

DOI: 10.1021/acs.jpcb.5b09213 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

overtone and combination bands. Third, it has been tested to the extent possible against higher-level quantum calculations for the water dimer and trimer and shown to produce energies within 15 cm−1 or less of the benchmark results.30 Finally, the comparison of the IR spectrum for liquid water at 300 K obtained with this method, using the WHBB potential and the current DMS, with experiment is very good, including the wellknown combination band at roughly 2200 cm−1.13

dipole moment using the DMS. The brackets represent an ensemble average of a NVT trajectory. For each temperature, 100 000 snapshots are used to obtain the average. The abovementioned ZPE-averaged monomer geometry is used for these calculations. 2.5. IR Spectrum. The IR spectrum is calculated using the local-monomer method.9 In brief, the Schroedinger equation for a given monomer configuration, [Tm̂ + Um(Q m) − Ei(m)]ϕi(m)(Q m) = 0

3. RESULTS AND DISCUSSION 3.1. Monomer Dipole Moment and Effective Monomer Charges. First, we present the distribution of monomer dipoles for the liquid simulation at 300 K in Figure 1 using the

(6)

is solved, where the T̂ m is kinetic operator and Um is the potential energy of monomer m with other monomers fixed. The Qm is the local-normal modes (lormal) of monomer m. These lormal modes are obtained from a lormal mode analysis, and there are nine modes for each monomer: three intramolecular ones and six intermolecular ones. For the current application, we use the q-TIP4P/F potential.3 The Schroedinger equation is solved using the VSCF/VCI method using the code MULTIMODE.26−28 The dipole transition moment is calculated by x ,y,z

R 0→f = m

∑ |⟨fm |Mα|0m⟩|2

(7)

α

0m and f m are the initial and the final states. Mα are components of the total dipole moment. In the calculation of the IR transition moment only the variation of the dipole moment with respect to local modes of monomer m needs to be considered, and this greatly simplifies the calculation of these matrix elements, which is done rigorously using numerical quadrature. As mentioned already, the published fifth order two-body DMS and a new, more efficient but somewhat less precise fourth-order DMS are used here. The IR spectrum is given by T

I (ν ) ∝

Figure 1. Distribution of monomer dipole magnitude of liquid water. See text for details.

two definitions above for this quantity. For visualization purposes the distributions are Gaussian-smoothed, as described above. The distribution from definition 1 is broad and peaked at 2.94 D. The average value is 2.94 D, which is roughly 1 D larger than the value of a gas phase water molecule. It is clear that the enhancement of monomer dipole is due to the twobody dipole terms. This result also agrees well with the value 2.9 D, inferred from X-ray diffraction experiments.15 The same experiment also suggests a charge transfer of 0.5 e along each OH bond. The monomer dipole distribution using the second definition is qualitatively similar to the one just discussed but shifted to lower values by around 16% and narrower. The average monomer dipole moment, 2.94 D, is also consistent with various previous theoretical predictions of 2.4−3.1 D using DFT, MP2, and empirical models.6,16,17 This Gaussian-like distribution agrees with previous results. Specifically, Sharma et al.31 obtained such a distribution based on ab initio (DFT) molecular dynamics calculations with 64 D2O molecules and localized Wannier functions for liquid water at 330 K. The peak in the distribution is at around 3.1 D. Shi et al. obtained such a distribution for liquid water at 303 K using their semiempirical approach to get the monomer dipole moment, with a peak at around 2.96 D. These results agrees well with the ab initio result shown here, using definition 1 for the monomer dipole moment. The emphasis here is not on the quantitative agreement, because different electronic structure theories would yield different results and the definition of this quantity is not unique across each study. Still, there is good consistency with these previous reports, suggesting again that truncating the

N

∑ ∑ ∑ ν0 → f (t ) R 0 → f (t ) Wf m

t = 0 m = 1 fm

m

m

g (ν−ν0 → f ) m

(8)

where the Wf m is a weighting factor defined as follow used for the LMon-4 approximation. This approximation has been applied previously,13,29 and a detailed description of this method has been given elsewhere.13 In brief, the approximation is to use six sets of 4-mode coupled calculations. In each set of calculations four modes are coupled. These are the three intramolecular modes of each monomer plus each of the six intermolecular modes (the frustrated rotations and translations) one-by-one. In this way, all nine lormal modes are considered. However, simply adding the eigenstates of all of the six sets of calculations would overcount the intramolecular states. To correct this, we use weighting factors of 1 for 6 intramolecular states and 1 for intermolecular states. More details on the weighting factor can also be found in our previous paper.13 Note that for IR spectrum calculation, the monomer sampling is from a 300 ps NVT trajectory, exactly the same sampling in the previous application on liquid water IR spectrum.13 This approach is clearly not exact; however, several points are worth making about it. First, and a trivial one, it is exact for the water monomer. And using the current one-body potential and dipole moment agreement with experiment is essentially perfect, in both position and intensity of the lines.13 Second, it does describe nonlinear aspects of the IR spectrum, i.e., D

DOI: 10.1021/acs.jpcb.5b09213 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

definition 2 to be roughly 20% lower, similar to the liquid water situation. There are 192 monomers in the ice model. We divide them into two groups: 105 monomers in the core region and 87 monomers on the surface. We show the monomer dipole distributions for the two groups of monomers. For the core monomers, the average monomer dipole is 3.54 D. This result is in apparent good agreement with the value of 3.56 D calculated (with an empirical scale factor) by Skinner and coworkers17 for ice Ih at 245 K. However, the large differences in temperature make a direct comparison problematic. Other theoretical studies predicted a range of values for ice Ih, from 2.3 to 3.7 D.31−36 The average monomer dipole for surface monomers is 3.02 D. To our best knowledge, this is the first time the surface monomer dipole is reported. The effective monomer charge distribution for ice Ih is shown in Figure 4.

many-body representation of dipole moment at the two-body level is realistic for liquid water. We analyzed the effective charges from the monomer dipole moments using eq 4, and the results are shown in Figure 2. By

Figure 2. Distribution of effective monomer charges of liquid water. See text for details.

use of the first definition, the average effective charge is 0.51 e on hydrogen atom and −1.02 e on oxygen atom. By use of the second definition, the average effective charge is 0.42 e on hydrogen atom and −0.84 e on oxygen atom. These results are qualitatively consistent with the experimentally suggested charge transfer level. Again, note that the effective charges in liquid water molecules are significantly different from the values of gas phase molecules, which are 0.33 e on hydrogen and −0.66 e on oxygen. This result reflects the change of local electronic environment from gas phase to liquid and is mathematically contributed from the two-body dipole terms. Also, note that in the q-TIP4P/F potential, the empirically parametrized hydrogen charge, 0.556 e, is also consistent with the present calculations. The monomer dipole distribution for ice Ih at 0 K is shown in Figure 3. To avoid redundancy in the figure, we only show the results from definition 1. One can expect the results from

Figure 4. Distribution of effective monomer charges of ice Ih. See text for details.

For core monomers, the average charge is 0.615 e on hydrogen and −1.22 e on oxygen. For surface monomers, the average charge is 0.504 e on hydrogen and −1.01 e on oxygen. The wider distribution of monomer dipole and monomer charges for surface monomers, compared to core monomers, indicates the higher heterogeneity on the surfaces. 3.2. Dielectric Constant. Since rigid monomers are used in the simulation, we begin by showing the sensitivity of the dielectric constant to the geometry of the rigid monomers. As mentioned in section 2, we use the ZPE-averaged geometry, and this results in a ϵ of 67 at 300 K. If simply using the gasphase equilibrium geometry, the ϵ would be 65. The ZPEaveraged result is closer to the experiment value of 78. The ZPE-averaged ϵ at the five temperatures, 280, 300, 320, 340, 360 K, are shown in Figure 5 along with the measured results.37 The dielectric constant shows a decreasing trend as temperature increases. This is very much expected, since the temperature and the volume, which expands as temperature increases, are in the denominator in eq 5. The dielectric constants are calculated using our DMS which contains both the one-body and two-body dipole terms. We have demonstrated in an earlier paper that the two-body dipole is essential to obtain correct IR intensity.13 It is also worthwhile to inspect the role of the two-body dipole in computing the dielectric constants. At 300 K the dielectric constant is 67 using the DMS(1b+2b). If we eliminate the two-body terms and use

Figure 3. Distribution of monomer dipole magnitude of ice Ih. See text for details. E

DOI: 10.1021/acs.jpcb.5b09213 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. Dielectric constant of liquid water at different temperatures.

only the one-body terms, the value is only 42, 37% smaller than the 1b+2b value. This demonstrates that the two-body dipole is very important for describing the dielectric constant or the total dipole of bulk water (or a large water cluster). Note that the two-body terms in our DMS are intrinsic two-body dipoles fit to ab initio data; therefore the electronic effects, such as polarization and charge transfer, are implicitly included in the fit. The result is consistent with the findings of Skinner and coworkers, that the original E3B model without polarization gives a dielectric constant of 54.2 at 303 K, but including polarization with a empirical linear model increases the value by 21.7.17 We do notice that the calculated values are about 10−15% lower than the experimental ones at each temperature. The following is several possible sources that may cause this difference. (1) Our many-body DMS is truncated at the twobody level. Previous tests7 show the 1b+2b dipole could be roughly 10% off the full ab initio dipole for water clusters up to 22-mer. (2) We assume vacuum boundary conditions when we evaluate the total dipole moment of the 500-mer cluster at each frame. The 500-mer may not be large enough to get a dipole that represents the situation in bulk water. The finite-size effect in periodic boundary conditions has been observed previously by Gereben and Pusztai.38 (3) Rigid monomers are assumed in the present simulations. (4) The dielectric constant also depends indirectly on the potential used in the MD simulation, although not significantly. We found that the value changes to 58 if using TIP4P potential with all other parameters same. Given that sensitivity to the potential, using accurate, flexible potentials, such as the ab initio WHBB water potential,7 would possibly predict a more accurate dielectric constant, albeit at a cost of significantly more computer time. Finally, we find that a 10 ns NVT trajectory is long enough for a converged dielectric constant. This agrees with Gereben and Pusztai, who claim at least 7 ns is needed for a converged result.38 3.3. IR Spectrum of Liquid Water. We show the IR spectrum of water at 300 K from LMon-4 calculations, using the q-TIP4P/F potential and three dipole moment surfaces, the fixed charge q-TIP4P/F one, the published ab initio, and also a lower order (fourth order) fit in Figure 6. The experimental spectrum at 298 K is also shown as a reference.39 The spectrum using the q-TIP4P/F fixed-charge DMS is shown in the top panel along with the spectrum reported previously using thermostated ring polymer molecular dynamics (TRPMD)29 using the q-TIP4P/F potential and dipole

Figure 6. IR spectra using q-TIP4P/F potential and various dipole moment surfaces. See text for details.

moment surface. First, note the good agreement between the LMon-4 and TRPMD spectra in both position and intensity of the major bands. The point of showing the comparison is to demonstrate the consistency of results using two very different methods (and different potentials used to obtain the thermal ensembles) over the large spectral range shown. (The agreement between the methods was also shown recently,29 although the LMon-4 method used in that study considered fewer intermolecular modes than the method used here.) There are some differences as well that are worth noting. The small peaks in the LMon spectrum at around 2350 and 3400 cm−1 are not seen in the TRPMD spectrum. The first one is due to combination band of the monomer bend plus an intermolecular libration. The second one is due to the overtone of the monomer bend. These are difficult to resolve in the TRPMD method. As seen, the calculated spectra are generally in good agreement with experiment in terms of the band positions. The q-TIP4P/F potential was optimized using RPMD to achieve this good agreement.3 Even so, the peak in the OHstretch band is roughly 60 cm−1 up-shifted from the experiment centered at 3413 cm−1. The broad band from 0 to 1000 cm−1 is the intermolecular band, characterizing the hindered translaF

DOI: 10.1021/acs.jpcb.5b09213 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

it uses configurations from classical MD simulations, currently based on efficient empirical water potentials like TIP4P and E3B. It is interesting to note that the approach of using a mixture of potentials and dipole moment surfaces in these simulations was advocated in a recent review article,40 and we clearly agree that this is a very useful strategy. The caveat to this approach is that the detailed results will certainly depend on the choice of the efficient potential and also on whether monomers are rigid or flexible. It is beyond the scope of this article to investigate this in detail. Some aspects of this were investigated in the recent paper by Rossi et al.,29 and further investigation is certainly warranted.

tional and rotational motions of water molecule. The band at roughly 1600 cm−1 is due to the intramolecular bending modes. It is interesting to note that the weak LMon combination band at 2350 cm−1, which is missing in the TRMPD spectrum mentioned above, is seen experimentally at around 2200 cm−1. Thus, the LMon peak is roughly 150 cm−1 upshifted relative to experiment but in “agreement” with the sum of the peak in the monomer bend (roughly 1690 cm−1) plus the peak in the libration band (roughly 700 cm −1 ). By contrast the experimental peak is roughly 150 cm−1 below the sum of the peak values seen experimentally. Is this an “error” in experiment? We doubt it. In fact, it is a deficiency in the qTIP4/F potential. The evidence for this is that when the WHHB PES is used in LMon calculations, the peak position of this combination band is only roughly 20 cm−1 up-shifted from experiment.13 The peak due to the bend overtone is not seen experimentally, and again this is a deficiency in the q-TIP4/F potential. The bend overtone does indeed contribute to the “OH-stretch” band; however, in previous calculations using the WHBB potential, this feature manifests itself in the breadth of the band and not as a distinct second peak.13 However, and this is the main point of this comparison, the intensities of the bands are not in good agreement with experiment. This is due to the limitation of fixed point-charge dipole models, which has been discussed previously.3,13 By contrast, the calculated LMon spectrum using the fifth order DMS shows significant improvement of the intensities and the agreement with the experimental spectrum is good. Because the DMS works very well along with the WHBB PES for the spectrum,13 it is not surprising that the ab inito DMS “corrects” the intensity from the q-TIP4P/F potential. Perhaps a little unexpected is that using the DMS also causes the band center of OH stretch shifts from 3475 cm−1 to 3451 cm−1. Given the fact that both spectra share the same vibrational density of states, this result suggests the dipole plays a role in the overall line shape of the spectrum by affecting the transition dipole moments, turning up and down “bright states”. The major difference between the ab initio DMS and the fixedcharge one is the increase in the dipole derivative in the Hbond OH-stretch mode, as demonstrated in detail previously.13 Finally, the spectrum using the more efficient fourth-order DMS, shown in the bottom panel, almost reproduces the one using the fifth-order fit. It is worthwhile to comment on the computational efficiency of the three above-mentioned dipoles. The simple, fixed point-charge q-TIP4P/F model is of course extremely efficient, a main reason for its popularity. The computational cost using the ab initio DMS combined with the q-TIP4P/F potentials is more but modest when used in LMon calculations. Using a single CPU, the intensity calculations take roughly 70 h with the fifth-order DMS and 11 h with the fourth-order DMS. Using q-TIP4P/F DMS results in roughly 3 h of CPU time. The solutions of the LMon Schroedinger equation for all monomers and all snapshots using the qTIP4P/F potential take roughly 60 h.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 404-727-6592. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank National Science Foundation (Grant CHE1463552) for financial support. We also thank Dr. Mariana Rossi for sending the TRMPD IR spectrum in digital format. We thank Yicun Ni and Prof. James Skinner for sending their modified version of GROMACS code with E3B potential.

■

REFERENCES

(1) 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. (2) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction models for water in relation to protein hydration. Intermol. Forces 1981, 14, 331. (3) Habershon, S.; Markland, T. E.; Manolopoulos, D. E. Competing quantum effects in the dynamics of a flexible water model. J. Chem. Phys. 2009, 131, 24501. (4) Wu, Y.; Tepper, H. L.; Voth, G. a. Flexible simple point-charge water model with improved liquid-state properties. J. Chem. Phys. 2006, 124, 024503. (5) Fanourgakis, G. S.; Xantheas, S. S. Development of transferable interaction potentials for water. V. Extension of the flexible, polarizable, Thole-type model potential (TTM3-F, v. 3.0) to describe the vibrational spectra of water clusters and liquid water. J. Chem. Phys. 2008, 128, 74506. (6) Wang, L.-P.; Head-Gordon, T.; Ponder, J. W.; Ren, P.; Chodera, J. D.; Eastman, P. K.; Martinez, T. J.; Pande, V. S. Systematic improvement of a classical molecular model of water. J. Phys. Chem. B 2013, 117, 9956−72. (7) Wang, Y.; Huang, X.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Flexible, ab initio potential, and dipole moment surfaces for water. I. Tests and applications for clusters up to the 22-mer. J. Chem. Phys. 2011, 134, 94509. (8) Lodi, L.; Tennyson, J.; Polyansky, O. L. A global, high accuracy ab initio dipole moment surface for the electronic ground state of the water molecule. J. Chem. Phys. 2011, 135, 034113. (9) Wang, Y.; Bowman, J. M. Ab initio potential and dipole moment surfaces for water. II. Local-monomer calculations of the infrared spectra of water clusters. J. Chem. Phys. 2011, 134, 154510. (10) Wang, Y.; Bowman, J. M. IR spectra of the water hexamer: Theory, with inclusion of the monomer bend overtone, and experiment are in agreement. J. Phys. Chem. Lett. 2013, 4, 1104−1108. (11) Liu, H.; Wang, Y.; Bowman, J. M. Local-monomer calculations of the intramolecular IR spectra of the cage and prism isomers of HOD(D2O)5 and HOD and D2O ice Ih. J. Phys. Chem. B 2014, 118, 14124−14131.

4. SUMMARY AND CONCLUSIONS The paper demonstrates the transferability of the recently reported, second-generation, ab inito DMS for water. It can be combined with model water potentials to obtain significantly more accurate IR spectrum of liquid water, with only modest extra computational efforts. It can be also applied in the electronic properties of bulk water, such as the monomer dipole distribution and dielectric constant. To obtain these properties, G

DOI: 10.1021/acs.jpcb.5b09213 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (12) Liu, H.; Wang, Y.; Bowman, J. M. Quantum calculations of intramolecular IR spectra of ice models using ab initio potential and dipole moment surfaces. J. Phys. Chem. Lett. 2012, 3, 3671−3676. (13) Liu, H.; Wang, Y.; Bowman, J. M. Quantum calculations of the IR spectrum of liquid water using ab initio and model potential and dipole moment surfaces and comparison with experiment. J. Chem. Phys. 2015, 142, 194502. (14) Medders, G. R.; Paesani, F. On the interplay of the potential energy and dipole moment surfaces in controlling the infrared activity of liquid water. J. Chem. Phys. 2015, 142, 212411. (15) Badyal, Y. S.; Saboungi, M. L.; Price, D. L.; Shastri, S. D.; Haeffner, D. R.; Soper, a. K. Electron distribution in water. J. Chem. Phys. 2000, 112, 9206−9208. (16) Gubskaya, A. V.; Kusalik, P. G. The total molecular dipole moment for liquid water. J. Chem. Phys. 2002, 117, 5290−5302. (17) Shi, L.; Ni, Y.; Drews, S. E. P.; Skinner, J. L. Dielectric constant and low-frequency infrared spectra for liquid water and ice Ih within the E3B model. J. Chem. Phys. 2014, 141, 084508. (18) Kemp, D. D.; Gordon, M. S. An interpretation of the enhancement of the water dipole moment due to the presence of other water molecules. J. Phys. Chem. A 2008, 112, 4885−4894. (19) Silvestrelli, P.; Parrinello, M. Water molecule dipole in the gas and in the liquid phase. Phys. Rev. Lett. 1999, 82, 3308−3311. (20) Tainter, C. J.; Pieniazek, P. a.; Lin, Y.-S.; Skinner, J. L. Robust three-body water simulation model. J. Chem. Phys. 2011, 134, 184501. (21) Kumar, R.; Skinner, J. L. Water simulation model with explicit three-molecule interactions. J. Phys. Chem. B 2008, 112, 8311−8318. (22) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511. (23) Hoover, W. G. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (24) Hayward, J. A.; Reimers, J. R. Unit cells for the simulation of hexagonal ice. J. Chem. Phys. 1997, 106, 1518. (25) Mas, E. M.; Szalewicz, K. Effects of monomer geometry and basis set saturation on computed depth of water dimer potential. J. Chem. Phys. 1996, 104, 7606−7614. (26) Carter, S.; Culik, S. J.; Bowman, J. M. Vibrational self-consistent field method for many-mode systems: A new approach and application to the vibrations of CO adsorbed on Cu(100). J. Chem. Phys. 1997, 107, 10458. (27) Carter, S.; Bowman, J. M.; Handy, N. C. Extensions and tests of “multimode”: a code to obtain accurate vibration/rotation energies of many-mode molecules. Theor. Chem. Acc. 1998, 100, 191−198. (28) Bowman, J. M.; Carter, S.; Huang, X. MULTIMODE: A code to calculate rovibrational energies of polyatomic molecules. Int. Rev. Phys. Chem. 2003, 22, 533−549. (29) Rossi, M.; Liu, H.; Paesani, F.; Bowman, J.; Ceriotti, M. Communication: On the consistency of approximate quantum dynamics simulation methods for vibrational spectra in the condensed phase. J. Chem. Phys. 2014, 141, 181101. (30) Wang, Y.; Bowman, J. M. Towards an ab initio flexible potential for water, and post-harmonic quantum vibrational analysis of water clusters. Chem. Phys. Lett. 2010, 491, 1−10. (31) Sharma, M.; Resta, R.; Car, R. Dipolar correlations and the dielectric permittivity of water. Phys. Rev. Lett. 2007, 98, 247401. (32) Rick, S. W. Simulations of ice and liquid water over a range of temperatures using the fluctuating charge model. J. Chem. Phys. 2001, 114, 2276−2283. (33) Rahman, A.; Stillinger, F. H. Proton distribution in ice and the Kirkwood correlation factor. J. Chem. Phys. 1972, 57, 4009−4017. (34) Heggie, M. I.; Latham, C. D.; Maynard, S. C.; Jones, R. Cooperative polarisation in ice Ih and the unusual strength of the hydrogen bond. Chem. Phys. Lett. 1996, 249, 485−490. (35) Gregory, J. K.; Clary, D. C.; Liu, K.; Brown, M. G.; Saykally, R. J. The water dipole moment in water clusters. Science 1997, 275, 814− 817. (36) Batista, E. R.; Xantheas, S. S.; Jónsson, H. Electric fields in ice and near water clusters. J. Chem. Phys. 2000, 112, 3285−3292.

(37) Fernández, D. P.; Mulev, Y.; Goodwin, A. R. H.; Sengers, J. M. H. L. A database for the static dielectric constant of water and steam. J. Phys. Chem. Ref. Data 1995, 24, 33−70. (38) Gereben, O.; Pusztai, L. On the accurate calculation of the dielectric constant from molecular dynamics simulations: The case of SPC/E and SWM4-DP water. Chem. Phys. Lett. 2011, 507, 80−83. (39) Bertie, J. E.; Lan, Z. Infrared intensities of liquids XX: The intensity of the OH stretching band of liquid water revisited, and the best current values of the optical constants of H2O(l) at 25 C between 15,000 and 1 cm−1. Appl. Spectrosc. 1996, 50, 1047−1057. (40) Vega, C. Water: one molecule, two surfaces, one mistake. Mol. Phys. 2015, 113, 1145−1163.

H

DOI: 10.1021/acs.jpcb.5b09213 J. Phys. Chem. B XXXX, XXX, XXX−XXX