Semiempirical Self-Consistent Polarization Description of Bulk Water

Mar 3, 2011 - In this work, we present the SCP-NDDO predictions for bulk water. ... Cartesian multipoles of ranks 0, 1, and 2 refer to the charge, dip...
0 downloads 0 Views 920KB Size
Download PDF
ARTICLE pubs.acs.org/JPCA

Semiempirical Self-Consistent Polarization Description of Bulk Water, the Liquid-Vapor Interface, and Cubic Ice Garold Murdachaew,† Christopher J. Mundy,* and Gregory K. Schenter† Chemical & Materials Sciences Division, Pacific Northwest National Laboratory, Richland, Washington 99352, United States

Teodoro Laino† Physical Chemistry Institute, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland, and IBM Zurich Research Laboratory, S€aumerstrasse 4, CH-8803 R€uschlikon, Switzerland

J€urg Hutter† Physical Chemistry Institute, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland ABSTRACT: We have applied an efficient electronic structure approach, the semiempirical selfconsistent polarization neglect of diatomic differential overlap (SCP-NDDO) method, previously parametrized to reproduce properties of water clusters by Chang, Schenter, and Garrett [J. Chem. Phys. 2008, 128, 164111] and now implemented in the CP2K package, to model ambient liquid water at 300 K (both the bulk and the liquid-vapor interface) and cubic ice at 15 and 250 K. The SCP-NDDO potential retains its transferability and good performance across the full range of conditions encountered in the clusters and the bulk phases of water. In particular, we obtain good results for the density, radial distribution functions, enthalpy of vaporization, self-diffusion coefficient, molecular dipole moment distribution, and hydrogen bond populations, in comparison to experimental measurements.

1. INTRODUCTION Despite large improvements in the algorithms of ab initio methods and advances in software and hardware, there will continue to exist system sizes and time durations for which the scaling of ab initio methods like CCSD(T), MP2, and even DFT remain prohibitive.1 Our own interest is in modeling the molecular condensed phase using molecular dynamics simulations, both with classical and quantum treatment of the nuclei. Because we would like to accurately describe many-body polarization effects and, eventually, the breaking and formation of chemical bonds, effective two-body classical potentials will not suffice, and we have to use a method that will faithfully describe the electronic structure during the course of the simulation. Most such direct dynamics or first-principles simulations2,3 have been performed using Kohn-Sham density functional theory (DFT),4,5 a method that is efficient enough to sample the configurational space, at least for some applications. Semiempirical (SE) methods have a long history in computational quantum chemistry. One class of these methods was originally developed to reduce the cost of the Hartree-Fock molecular orbital method and, through parametrization, to include some of the missing correlation effects. The modern version of this approximation is based on neglect of diatomic differential overlap (NDDO).6,7 Some widely used parametrizations have been the MNDO,8 AM1,9 PM3,10,11 and PM612 approaches. The other class of SE methods is the tight-binding approximation to DFT, for example, the self-consistent charge r 2011 American Chemical Society

density functional tight-binding (SCC-DFTB)13 method. SE methods are less costly than DFT. Therefore, they may enable more extensive sampling during a simulation. Both NDDObased and DFTB-based SE methods have been an effective and efficient approach for chemically bonded complexes. Some deficiencies of DFT are the lack of the dispersion interaction and an often poor description of weakly bound systems.14 Thus, new exchange-correlation functionals15,16 and empirical dispersion corrections17,18 have been developed for DFT. Similar corrections have also been applied to the NDDObased19,20 and to the DFTB-based21,22 SE methods. A general conclusion is that no purely ab initio method exists for an on-thefly electronic structure description that is both efficient and accurate. The two currently feasible approaches, DFT and SE (either NDDO or DFTB), both require some degree of tuning to accurate experimental or ab initio data. Here, we will focus on the NDDO-based SE methods. There have been many attempts to improve the NDDO methods’ description of the hydrogen bond. A pairwise core repulsion function was added to MNDO in the MNDO/H,23 MNDO/ GH,24 MNDO/M,25,26 and I-MNDO27 methods. After MNDO, Special Issue: Victoria Buch Memorial Received: November 2, 2010 Revised: January 25, 2011 Published: March 03, 2011 6046

dx.doi.org/10.1021/jp110481m | J. Phys. Chem. A 2011, 115, 6046–6053

The Journal of Physical Chemistry A

ARTICLE

the newer AM1,9 PM3,10 and PM612 parametrizations were introduced. To improve the description of hydrogen bonds, additional empirical functions where added in the (AM1,PM3)PIF,19 PM3-MAIS,20 and PDDG-(MNDO,PM3)28 methods. A different definition of the electron integrals was used in the SAM1 method.29 Additional basis functions of p symmetry were added to H atoms in the SINDO/130 method. Orthogonalization corrections were introduced in the OMn31,32 methods. As mentioned, post-SCF dispersion corrections1,33 have also been added. The NDDO-based SE methods have been used to model the condensed phase (see, e.g., refs 34-38), as have the DFTBbased SE methods. For simulations of ambient water, see, for example, refs 21, 22, and 38. In ref 39, some of us introduced a different correction approach. The self-consistent polarization (SCP) approach uses an auxiliary charge density (modeled by one or at most a few Gaussian centers placed on atomic sites) which supplements the NDDO or DFT charge densities and is used to better describe the intermolecular interaction. SCP-NDDO was developed first and used to correct the deficiencies of the PM3 description of water clusters39 (geometries, binding energies, and harmonic frequencies). Later, the same SCP approach was used to correct the deficiencies of the BLYP DFT functional in describing argon clusters40 and water clusters.41 The SCP-NDDO and SCP-DFT approaches were parametrized on properties of clusters and used to predict other cluster properties and, in a successful test of transferability, bulk properties, in good agreement with accurate theoretical and experimental benchmarks (SCP-DFT for liquid and solid argon,40 radial distribution functions, or RDFs, and cohesion energy curves; SCP-DFT for liquid ambient water,41 RDFs, and enthalpy of vaporization). In this work, we present the SCP-NDDO predictions for bulk water. Water in the liquid, interfacial, and solid states has been of interest to many fields, and ambient liquid water, in particular, has been studied computationally for decades with a variety of potentials.42 An efficient and accurate description which includes the electronic state would be of interest to many fields including biology43 and catalysis and atmospheric science.44 In this work, we present our SCP-NDDO results for ambient bulk and interfacial water and also for cubic ice, as obtained from SE molecular dynamics simulations.

2. IMPLEMENTATION OF NDDO AND SCP-NDDO FOR EXTENDED SYSTEMS 2.1. Overview of SCP-NDDO Theory. The SCP-NDDO method was fully described in ref 39. Some essential points are summarized here. The SCP-NDDO total energy is given by

E ¼ Eself þ Ecoul þ Eex þ Ecore þ Edisp

ð1Þ

where the respective terms’ subscripts refer to the self, Coulomb, resonance þ exchange, core, and dispersion contributions. The self and Coulomb terms have SCP contributions which differentiate the SCP-NDDO method from the parent MNDO and PM3 methods. Many of the parameters have been refitted39 and thus differ slightly from their original values.10 In addition, the screening function (and thus the integrals) in the terms differ from those in the parent methods. Diagonalization of the associated Fock matrix yields the self-consistent problem, which is solved variationally for the electron and polarization density

matrices. The dispersion term is an additional post-SCF (i.e., postprocessed) term. In SCP-NDDO, each product of two atomic orbitals (AOs) at a single atomic center occurring in a general two-center electron repulsion integral is approximated by a Cartesian point multipole. Thus, SCP-NDDO differs from standard NDDO methods in that those methods further write the multipoles as distributed charges around the center. We note also that the NDDO (e.g., MNDO, AM1, PM3, and PM6) methods are programmed in CP2K with atomic point multipoles rather than distributed charges. Advantages of using point multipoles at atomic sites rather than charges only, at atomic and at additional off-atomic sites, is ease of implementation and increased efficiency. Because there are, at most, AOs of p symmetry on atoms, the maximum multipole needed is the quadrupole. Thus, a general two-electron integral can be approximated as ð0Þ

ð1Þ

ð2Þ

ðμν; AAjλσ; BBÞ  ðM μν;A þ M μν;A 3 rA þ M μν;A : rA rA Þ ð0Þ

ð1Þ

ð2Þ

ðM λσ;B þ M λσ;B 3 rB þ M λσ;B : rB rB Þ sðRAB Þ  RAB

ð2Þ

where the subsequent Cartesian multipoles of ranks 0, 1, and 2 refer to the charge, dipole, and quadrupole, respectively, and we have introduced the screening function present in NDDO theories. Each multipole is analytically precoded, and the set of multipoles need to be calculated just once for each system (i.e., given the atomic kinds and the corresponding AO basis sets) using the standard definitions. See ref 39 for additional details. 2.2. Multipolar Ewald Summation Approach and Screening Function. Implementations of NDDO methods for periodic systems have appeared in the literature.34-36 References 37 and 38 present examples of calculations on ambient water using various approaches to treat the extended system. Here, we describe the implementation of SCP-NDDO within CP2K,45 which was part of a broader effort to code NDDO-based SE methods. There are two major complications with implementing NDDO methods for periodic systems; the screening function may cause problems for certain systems,35,36 and if multipole expansions are used instead of point charges, there is the need to code multipolar Ewald46,47 summations. In SCP-NDDO, the first problem was addressed by use of a more appropriate screening function. For all NDDO methods, the second problem was resolved by efficiently coding multipolar Ewald summations and consistent gradients up to the quadrupolar term using the derivations in refs 48 and 49. In standard SE NDDO-based methods like MNDO, AM1, PM3, and PM6, the Klopman-Dewar-Sabelli-Ohno (KDSO)50-54 screening function is used, and sKDSO ðRAB Þ 1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RAB 2 RAB þ ðFA þ FB Þ2

ð3Þ

where FA and FB are screening parameters for atomic kinds A and B, respectively. In ref 39, we elected to develop SCP-NDDO with the more physically motivated Slater screening function, 6047

dx.doi.org/10.1021/jp110481m |J. Phys. Chem. A 2011, 115, 6046–6053

The Journal of Physical Chemistry A and thus sSlater ðRAB Þ ¼ RAB

1 - eðσ

A

þ σB ÞRAB

ARTICLE

  1 1 þ ðσ A þ σB ÞRAB 2 RAB

ð4Þ

where σA and σB are screening parameters for atomic kinds A and B, respectively. For both screening functions, we have used Ewald techniques. The Ewald separation of the operator into a reciprocal and a rapidly converging real space sum is facilitated by first subtracting and adding the Coulomb operator     sðRAB Þ sðRAB Þ 1 1 ¼ þ ð5Þ RAB RAB RAB RAB The second term in eq 5 is just the Coulomb operator, and in eq 2, it can be handled using the multipolar Ewald summation procedure. For the Slater screening, it is obvious that the first term in eq 5 is explicitly short-range (it decays exponentially) and thus can be handled by a direct space sum with a small cutoff (in practice, it is calculated together with the short-range/direct space part of the Ewald term). If using the KDSO screening, the first term also contributes a leading-order dipolar term decaying as (1/R3AB), something which can be readily shown from a bipolar expansion. A dipolar Ewald summation could be applied to this as well, but in practice, we found that this was not necessary (the error introduced in neglecting the dipolar Ewald summation of this term was negligible). 2.3. Charge Constraint for the Condensed Phase. The SCP-NDDO parameters used were identical to the ones fitted39 to properties of water clusters. However, in some preliminary simulations on cubic ice at 100 K and on ambient water, we observed that a few molecules were undergoing unphysical dissociation. In the cases where the dissociation was occurring, it was observed that the molecules were becoming overly polarized, as reflected in the Mulliken charge QA of the Ath atom, where QA = ∑μ(pμμ;A - PTμμ;AA), the total SCP plus NDDO electronic charge on atom A. The only changes made to the SCPNDDO water Hamiltonian obtained in the earlier work39 was to add a penalty function to prevent an atom’s charge from going outside of specific bounds. The piecewise harmonic charge constraint penalty function was defined as ΔE

¼ ¼

∑A 12kA ðQA - QA;0 Þ2 ∑A 12kA ðQA - QA;1 Þ2

QA < QA;0

ð6Þ

QA > QA;1

To ensure self-consistency in the SCF minimization, it was necessary to compute the contribution ΔFμμ;AA to the Fock matrix Fμν;AB, which was ΔFμμ;AA

¼ kA ðQA - QA;0 Þ ¼ kA ðQA - QA;1 Þ

QA < QA;0 QA > QA;1

ð7Þ

The constants kA and the lower and upper charge limits QA;0 and QA;1 (valence only, by convention taken to be positive) needed to be defined for each atom type. After some tests on the ambient liquid, we chose kO = kH = 1.0, QO;0 = 6.1, QO;1 = 6.6, QH;0 = 0.72, and QH;1 = 0.92. The charge constraint was implemented in a local version of the CP2K code. We confirmed that the CP2K implementation of SCP-NDDO gave identical results for clusters (nonperiodic calculation) with or without the penalty function

and was in agreement with the earlier work.39 We note that the charge constraint is a temporary correction. Better short-range (i.e., intramolecular) screening in the future should obviate the need for such a constraint in the condensed phase. 2.4. Cost of Various Simulation Methods. In using CP2K to perform NVT simulations of ambient water with various methods, we confirmed that the NDDO-based SE methods (PM3, PM6, and SCP-NDDO) are more than an order of magnitude cheaper than DFT. Just like DFT, the SE methods scale linearly with the number of processors. Work is underway to make the implementation even more efficient.

3. COMPUTATIONAL DETAILS The CP2K package45 was used for all simulations. The SE PM3, PM6, and SCP-NDDO calculations were performed using the semiempirical module as implemented in CP2K. The BornOppenheimer molecular dynamics simulations used a time step of 0.5 fs. The duration of the water and ice canonical (NVT) simulations was 40 ps (10 ps for equilibration with the target potential and 30 ps for production of averages). To check convergence in time, we confirmed that using the first 10 ps of production gave essentially the same average properties as those based on the full 30 ps production run. To better quantify the statistical errors, we followed the block averaging approach of ref 55; the 30 ps simulation was subdivided into three equal blocks, and block averages were obtained. The standard deviation of the block averages from the mean over the full simulation is presented as our statistical uncertainty. The ambient liquid water NVT simulations were performed with 64 water molecules in a periodic cubic box of length L = 12.4138 Å, resulting in a density of F = 1.0 g/cm3. The wave function convergence criterion was 10-7 au. The temperature was held constant at 300 K using the method of Nose and Hoover,56-59 with chain thermostats applied to all degrees of freedom (so-called “massive” thermostatting60,61). The starting configuration was a well-equilibrated configuration obtained from previous work in our group. Whereas static properties (e.g., radial distribution functions, internal energy/enthalpy of vaporization, molecular dipole distribution, and hydrogen bond populations) at a given temperature can be obtained using the canonical ensemble, the dynamical properties (e.g., self-diffusion coefficient and time correlation functions and the related infrared and neutron spectra) are best obtained using the microcanonical (NVE) ensemble, where energy is strictly conserved. Thus, we performed NVE simulations in the following manner. For the SCP-NDDO ambient liquid, 21 starting configurations were selected, and for each configuration, 10 ps long simulations were performed in a manner identical to the NVT parent simulation, with the exception that the thermostats were disabled. The 21 starting points were the 5th through the 25th ps points of the parent NVT simulation. This resulted in having a total of 210 ps of an NVE trajectory. As could have been expected, our examination of the SCP-NDDO data showed that the most slowly convergent property was the self-diffusion coefficient. Thus, for the systems and potentials (e.g., PM3 and PM6) for which a less precise estimate of this property was sufficient, we performed just five 10 ps NVE simulations whose starting points were the 5th, 9th, 13th, 17th, and 21st ps points of the parent NVT trajectory. The multiple starting points allowed for block averaging and determination of errors, as described above. Although there was some 6048

dx.doi.org/10.1021/jp110481m |J. Phys. Chem. A 2011, 115, 6046–6053

The Journal of Physical Chemistry A natural temperature variation among the NVE trajectories, the average temperature over the total NVE trajectory for each potential was close to the target temperature. The SCP-NDDO water slab NVT calculations needed to model the ambient liquid-vapor interface were performed in a similar fashion to the bulk water NVT calculations, with the exceptions that 216 water molecules were used and the dimensions of the rectangular box were 15 Å (x)  15 Å (y)  71.44 Å (z). Such a setup had been found62,63 to be the minimal system size needed to ensure the presence of a stable interfacial system with a bulk-like water region with a width of about 26 Å and a vacuum/vapor region of about 23 Å on either side of the periodic slab. As in those previous calculations,62,63 a constraint was used to ensure that the center of mass of the slab remained at the origin of the coordinate system. This served to simplify the subsequent analysis. The duration of the SCP-NDDO run was 14.8 ps (4 ps of equilibration and 10.8 ps of production). The starting configuration was a well-equilibrated configuration from our previous work63 using DFT. The water molecule was modeled as D2O rather than H2O in the slab simulations. Within the classical nuclei approximation, such a change of mass has no effect on static properties. SCP-NDDO cubic ice Ic NVT simulations were performed in a similar manner to the bulk liquid simulations, except that the wave function convergence criterion was tightened to 10-8 au for better energy conservation. We elected to model cubic ice Ic (rather than, e.g., hexagonal ice Ih) because this form permits the use of a small primitive unit cell which is ideal for test calculations. After performing some initial calculations using a primitive unit cell of side length L = 6.358 Å containing 8 water molecules, we doubled the cell in all dimensions to obtain a cubic cell of side length L = 12.716 Å containing 64 water molecules. The use of a larger unit cell allows a better sampling of the Brillouin zone in our Γ-point calculation. The density of our cubic ice simulations was F = 0.9309 g/cm3. The temperatures used were 250 and 15 K. The starting configurations were first geometry-optimized at 0 K. This was followed by a short NVE simulation at twice the target temperature. Finally, the NVT simulation at the target temperature was performed. For the bulk liquid and ice SE simulations, electrostatic interactions were treated using Ewald summation,48,49 with the real space sum truncated at L/2 (about 6.4 Å). For the slab calculations, the real space summation was truncated at 7.5 Å. For all simulations and potentials, the summation in reciprocal space included G vectors up to Gmax = 25 and an Ewald real/reciprocal space summation separation parameter of R = 0.5 Å-1.

4. RESULTS 4.1. Dipole moment distribution in the liquid. Table 1 summarizes our results for the average dipole moment of the monomer in the liquid. Gas-phase values are also shown, as are comparisons to literature results. For the NDDO-based methods, the monomer’s dipole moment was calculated as the sum of the atomic dipole and charge contributions readily extractable from our NDDO implementation. The benchmark values of the gas-phase (experimental) and liquid (theoretical) water molecular dipole moment are 1.85 D64 and about 2.6565 to 3 D,66 with an experimental study producing a value of 2.9 D.67 Figure 1 shows the molecular dipole moment distribution in ambient liquid water as predicted from the PM3, PM6, and SCP-NDDO simulations. In Table 1, the average

ARTICLE

Table 1. Dipole Momenta μ [D] Gas Phase expt

1.85b

PM3

1.74c

PM6

2.07

SCP-NDDO

2.16c Ambient Liquid

expt/thry PM3

2.9d/2.65e; 3f 1.9

PM6

2.3

SCP-NDDO

2.8

a

Presented is the dipole moment of the isolated molecule and its average value in the ambient liquid obtained from the NVT simulations and comparisons to the literature. For the ambient liquid, T = 300 K unless otherwise indicated. b Reference 64. c Reference 39. d Reference 67 (T = 298 K). e Reference 65 (T = 298 K). f Reference 66 (T = 298 K).

Figure 1. Molecular dipole moment distribution in ambient liquid water.

values are presented. The average PM3 value in the liquid (1.9 D) is essentially the same as that for the isolated PM3 molecule (1.74 D), whereas, the SCP-NDDO model gives a good result; the dipole moment increases from the isolated molecule value of 2.16 to 2.8 D in the liquid. The dipole moment predicted by PM6 is intermediate between PM3 and SCP-NDDO, 2.07 D for the isolated molecule and 2.3 D for the liquid. Thus, the auxiliary polarization density present in the SCP-NDDO method corrects for the lack of polarization in the minimum basis set NDDO methods. 4.2. Radial Distribution Functions in the Liquid and in Cubic Ice. The PM3, PM6, and SCP-NDDO RDFs gOO(r), gOH(r), and gHH(r) are shown in Figure 2. They are compared to an experimental benchmark, the neutron data of Soper.68 In Table 2, we summarize our liquid-state simulations by first presenting the first maximum, first minimum, and second maximum of gOO(r), and we compare to literature results. By integrating gOO(r) from the origin to the first minimum, the coordination number (CN) is obtained. First, consider gOO(r) in Figure 2 and in Table 2, in particular, the results given by the uncorrected SE method, PM3. While PM3 has the first maximum in the correct general location, the 6049

dx.doi.org/10.1021/jp110481m |J. Phys. Chem. A 2011, 115, 6046–6053

The Journal of Physical Chemistry A

ARTICLE

Figure 2. Radial distribution functions of ambient liquid water compared to the experimental data of Soper.68

Table 2. Results of the NVT and NVE Simulations of Ambient Water and Comparisons to the Literaturea

expt

max1 (rmax1 OO ,gOO )

min1 (rmin1 OO ,gOO )

max2 (rmax2 OO ,gOO )

([Å],[unitless])

([Å],[unitless])

([Å],[unitless])

(2.73,2.75)b

(3.36,0.78)b

(4.52,1.16)b

h

h

h

(2.73,2.83)

(3.41,0.79)

(4.44,1.13)

CN

H bonds

U [kcal/

ΔHvapb

mol]

[kcal/mol]

3.58d; 3.5e -9.9 ( 0.3f

4.5b,c

TNVE [K]

10.5 ( 0.3f

D [Å2/ps] 0.23g

h

4.7

PM3

(2.70,3.33 ( 0.02) (2.97,0.86 ( 0.03) (3.33,1.24 ( 0.03) 2.75 ( 0.04 2.0

-4.00 ( 0.02

4.60 ( 0.02

303 ( 4 0.78 ( 0.08

PM6

(2.64,2.32 ( 0.02) (3.88,0.89 ( 0.03) (5.03,1.10 ( 0.03) 7.79 ( 0.04 1.5

-5.66 ( 0.02

6.26 ( 0.02

301 ( 4 1.28 ( 0.02

SCP-NDDO

(2.74,2.95 ( 0.07) (3.28,0.69 ( 0.03) (4.40,1.27 ( 0.03) 4.17 ( 0.04 3.3

-10.05 ( 0.05 10.65 ( 0.05

299 ( 3 0.16 ( 0.01 (0.21i)

a

Presented are the oxygen-oxygen RDF main features, the coordination number, the average number of hydrogen bonds per molecule, the internal energy and enthalpy, and also the self-diffusion coefficient. T = 300 K unless otherwise indicated. Error bars are given where available. The error on the position of the oxygen-oxygen RDF main features is estimated to be (0.02 Å. See the text for additional details. b Reference 68 (T = 298 K). c Taken from ref 70. d Reference 76 (T = 298 K). e Reference 77 (T = 298 K). f Reference 78 (T = 298 K). g Reference 71 (T = 298 K). h References 79 and 80. i With finite-system size correction.

subsequent features (first minimum and second maximum) are shifted inward, and overall, PM3 fails for predicting the structure. PM6 has features shifted dramatically outward and thus also fails in predicting the correct structure. The SCP-NDDO method, on the other hand, behaves reasonably and is close to the experimental curve. SCP-NDDO also gives good gOH(r) and gHH(r) curves. The SCP-NDDO gOO(r) is slightly overstructured compared to the experimental curve. Quantum nuclear motion corrections, based on differences between classical and quantum simulations, are known to reduce the overstructuring of gOO(r) curves obtained with classical treatment of the nuclei (see, e.g., refs 69 and 70). Thus, the quantum corrections would bring the SCPNDDO gOO(r) features into even better agreement with the experimental benchmark. Figure 3 shows the effect of phase and temperature on the SCP-NDDO RDFs. As expected, more structure is visible in the transition from liquid water at 300 K to cubic ice, first at 250 K and then at 15 K. 4.3. Enthalpy of Vaporization in the Liquid. The enthalpy of vaporization ΔHvap = -U þ RT (where R is the gas constant and

U is the internal energy) is presented in Table 2. The SCPNDDO result of 10.66 kcal/mol is in excellent agreement with the experimental value of 10.5 ( 0.3 kcal/mol. The PM3 and PM6 results are too small, reflecting the lack of sufficient polarization and dispersion interactions. 4.4. Self-Diffusion and Dynamics in the Liquid. For each of the 21 SCP-NDDO and 5 PM3 and PM6 NVE trajectories, the self-diffusion coefficient D was obtained from the Einstein relation 6Dt ¼ ÆjrðtÞ - rð0Þj2 æ

ð8Þ

where the term on the right-hand side is the mean-square displacement (MSD). The notation Æ 3 3 3 æ refers to averaging over all relevant quantities, in this case, the displacements of all oxygen atoms and multiple starting points separated from one another by 100 fs. This separation was chosen because this is approximately the velocity autocorrelation time (i.e., the time for the oxygen-oxygen velocity autocorrelation function of the various potentials to decay to 0). D was determined by fitting a straight line to the MSD at intermediate times (as in ref 55 to the 6050

dx.doi.org/10.1021/jp110481m |J. Phys. Chem. A 2011, 115, 6046–6053

The Journal of Physical Chemistry A

ARTICLE

Figure 3. Effect of phase and temperature on the SCP-NDDO RDFs.

Figure 4. Mean-square displacement as a function of time.

2-7 ps portion of each 10 ps trajectory) to avoid contamination from the ballistic regime at smaller times (in practice, the region with t j 300 fs) and the lower statistics region at the end of the trajectory. Averaging was then performed, and the average MSDs thus obtained from the NVE PM3, PM6, and SCP-NDDO simulations are shown in Figure 4, while in Table 2, we present the mean D values, with the statistical uncertainty shown being the standard deviation of the mean value over the NVE trajectories. The mean temperatures and their uncertainties are also presented because there was some small drift in the NVE simulations. The PM3 and PM6 diffusion coefficients are much too large, while the SCP-NDDO result of 0.16 Å2/ps is in reasonable agreement with the experimental value71 of 0.23 Å2/ps. It is known that the use of a small system size results in a smaller diffusion coefficient. Thus, calculations on systems of 64 and 512 molecules were performed using the SPC/E potential.72

Figure 5. Density, molecular dipole, and hydrogen bond population profiles as obtained from the slab calculation using the SCP-NDDO potential. The thick black line is the fit to the density profile.

The density was 1.0 g/cm3. While negligible differences were observed in most static properties (e.g., RDFs and cohesion energies), the SPC/E diffusion coefficient increased from 0.26 to 0.28 Å2/ps. Extrapolating D to infinite system size using eq 12 in ref 73, a value of 0.31 Å2/ps was obtained. Note that this is in good agreement with the SPC/E literature value74 of 0.249 Å2/ ps. Adding the estimated finite-size correction, the SCP-NDDO diffusion coefficient becomes 0.21 Å2/ps, and the agreement with experiment is improved. Note that this approximation relies on the assumption that the SCP-NDDO viscosity is close to that of SPC/E. 4.5. Density, Dipole, and Hydrogen Bond Populations at the Interface and in the Liquid. Figure 5 shows how the SCPNDDO density, dipole moment, and hydrogen bond populations behave as one transitions from the interior to the outer 6051

dx.doi.org/10.1021/jp110481m |J. Phys. Chem. A 2011, 115, 6046–6053

The Journal of Physical Chemistry A

ARTICLE

Table 3. Results of the NVT Simulations of the Ambient Liquid-Vapor Interface and Comparison to Experimenta Fl [g/cm3] expt

0.995

SCP-NDDO

0.993

zGDS [Å]

δ [Å]

14.43

1.010

b

a

Presented are the liquid density and the position and thickness of the Gibbs dividing surface. T = 300 K unless otherwise indicated. b Reference 68 (T = 298 K).

Figure 6. Like Figure 5 with SCP-NDDO hydrogen bond populations, but slab results are averaged over two regions; bulk results are also shown for comparison.

region of the slab. The data were obtained using a bin width of 1 Å and by averaging over both surfaces of the slab. As in refs 62, 63, and 75, the density profile was fitted using   1 1 z - zGDS FðzÞ ¼ ðFl þ Fv Þ - ðFl - Fv Þ tanh ð9Þ 2 2 δ where Fl and Fv are the liquid and vapor densities, respectively (the latter is taken to be 0); zGDS is the position of the Gibbs dividing surface (taken to be the point where the density is half of that of the liquid density), and δ is the interfacial width (related to the 10-90% thickness t by t = 2.197δ). The results of the fit are shown in Table 3. The SCP-NDDO liquid density of 0.993 g/ cm3 is in agreement with the experimental68 value of 0.995 g/ cm3. Note that such close agreement is likely to be at least partially accidental. Figure 6 shows similar data but subdivides the slab into two regions and also includes a comparison of the slab interior to the results from the bulk calculation. As in ref 63, the division of the slab into interior and surface regions was done by defining as interior molecules all of those molecules with oxygen positions z less than zGDS - 2δ. Our results are in good agreement with those obtained previously62,63 using a much more expensive DFT (BLYP) potential. Note that in the bulk or interior, there is mainly tetrahedral organization (the DDAA type dominates), whereas in the outer or vapor region, there is more ring-like organization (the DA type dominates). Thus, the water molecules in the interface region generally lie in the interfacial plane.

5. CONCLUSIONS The SCP-NDDO SE method previously parametrized39 for water clusters had been found to perform much better than PM3 and PM6. After implementing the NDDO-based SE methods in the CP2K package through an efficient coding of the multipolar Ewald approach, we performed simulations of bulk water with these methods. Just as for the clusters, the PM3 and PM6 methods performed poorly for the heat of vaporization, dipole moment, and RDFs in the ambient liquid, while the SCP-NDDO method yielded good agreement with experiment. Particularly in the case of RDFs of liquid water, this is a noteworthy result because in the past, SE methods (both NDDO- and DFTB-based,

even with empirical corrections) have struggled to even qualitatively reproduce experimental RDFs. Our simulations on cubic ice also yielded reasonable results. Simulations on the liquidvapor interface have yielded density, dipole moment, and hydrogen bond population profiles in good agreement with experiment and with earlier work using a more expensive DFT potential. We believe that in our SCP-NDDO implementation, the SE approach (either NDDO- or DFTB-based), which in the past failed for describing water in the bulk, now holds promise as a method that can accurately and efficiently be used in first-principles simulations of the condensed state.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Phone: þ1 509-375-2404. Fax: þ1-509-375-2644. Notes †

E-mail: [email protected] (G.M.); greg.schenter@ pnl.gov (G.K.S.); [email protected] (T.L.); [email protected] (J.H.).

’ ACKNOWLEDGMENT This work was performed under the auspices of the Division of Chemical Sciences, Geosciences, and Biosciences, Office of Basic Energy Sciences, U.S. Department of Energy, under Contract No. DE-AC06-76RLO 1830 with Battelle Memorial Institute, which operates the Pacific Northwest National Laboratory, a multiprogram national laboratory. This research was performed in part using the computational resources in the National Energy Research Supercomputing Center (NERSC) at Lawrence Livermore National Laboratory. We also used the NWICE computational resource at PNNL. G.M. would like to acknowledge discussions with Shawn Kathmann. ’ REFERENCES (1) Clark, T. J. Mol. Struct.: THEOCHEM 2000, 530, 1. (2) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471. (3) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge: New York, 2009. (4) Hohenberg, P.; Kohn, W. Phys. Rev. B 1964, 136, 864. (5) Kohn, W.; Sham, L. J. Phys. Rev. A 1965, 140, 1133. (6) Pople, J. A.; Santry, D. P.; Segal, G. A. J. Chem. Phys. 1965, 43, S129. (7) Sustmann, R.; Williams, J. E.; Dewar, M. J. S.; Allen, L. C.; Schleyer, P. v. R. J. Am. Chem. Soc. 1969, 91, 5350. (8) Dewar, M. J. S.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 4899. (9) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (10) Stewart, J. J. P. J. Comput. Chem. 1989, 10, 209. (11) Stewart, J. J. P. J. Mol. Model. 2004, 10, 155. 6052

dx.doi.org/10.1021/jp110481m |J. Phys. Chem. A 2011, 115, 6046–6053

The Journal of Physical Chemistry A (12) Stewart, J. J. P. J. Mol. Model. 2007, 13, 1173. (13) Elstner, M.; Porezag, D.; Jungnickel, G.; Elstner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58, 7260. (14) Wu, X.; Vargas, M. C.; Nayak, S.; Lotrich, V.; Scoles, G. J. Chem. Phys. 2001, 115, 8748. (15) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215. (16) Pernal, K.; Podeszwa, R.; Patkowski, K.; Szalewicz, K. Phys. Rev. Lett. 2009, 103, 263201. (17) Grimme, S. J. Comput. Chem. 2006, 27, 1787. (18) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (19) Bernal-Uruchurtu, M. I.; Martins-Costa, M. T. C.; Millot, C.; Ruiz-Lopez, M. F. J. Comput. Chem. 2000, 21, 572. (20) Bernal-Uruchurtu, M. I.; Ruiz-Lopez, M. F. Chem. Phys. Lett. 2000, 330, 118. (21) Hu, H.; Lu, Z.; Elstner, M.; Hermans, J.; Yang, W. J. Phys. Chem. A 2007, 111, 5685. (22) Maupin, C. M.; Arad, B.; Voth, G. A. J. Phys. Chem. B 2010, 114, 6922. (23) Burstein, K. Y.; Isaev, A. N. Theor. Chim. Acta 1984, 64, 397. (24) Goldblum, A. J. Comput. Chem. 1987, 8, 835. (25) Voityuk, A. A.; Bliznyuk, A. A. Theor. Chim. Acta 1987, 72, 223. (26) Bliznyuk, A. A.; Voityuk, A. A. J. Mol. Struct.: THEOCHEM 1988, 41, 343. (27) Salk, S. H. S.; Chen, T. S.; Hagen, D. E.; Lutrus, C. K. Theor. Chim. Acta 1986, 70, 3. (28) Repasky, M. P.; Chandrasekhar, J.; Jorgensen, W. L. J. Comput. Chem. 2002, 23, 1601. (29) Dewar, M. J. S.; Jie, C. X.; Yu, J. G. Tetrahedron 1993, 49, 5003. (30) Jug, K.; Geudtner, G. J. Comput. Chem. 1993, 14, 639. (31) Kolb, M.; Thiel, W. J. Comput. Chem. 1993, 14, 775. (32) Weber, W.; Thiel, W. Theor. Chem. Acc. 2000, 103, 495. (33) Winget, P.; Selcuki, C.; Horn, A. H. C.; Martin, B.; Clark, T. Theor. Chem. Acc. 2003, 110, 254. (34) Gale, J. D. Faraday Discuss. 1997, 106, 219. (35) Stewart, J. J. P. J. Mol. Struct. 2000, 556, 59. (36) Stewart, J. J. P. J. Mol. Model. 2008, 14, 499. (37) Monard, G.; Bernal-Uruchurtu, M. I.; van der Vaart, A.; Merz, K. M., Jr.; Ruiz-Lopez, M. F. J. Phys. Chem. A 2005, 109, 3425. (38) Geerke, D. P.; Thiel, S.; Thiel, W.; van Gunsteren, W. F. Phys. Chem. Chem. Phys. 2008, 10, 297. (39) Chang, D. T.; Schenter, G. K.; Garrett, B. C. J. Chem. Phys. 2008, 128, 164111. (40) Maerzke, K. A.; Murdachaew, G.; Mundy, C. J.; Schenter, G. K.; Siepmann, J. I. J. Phys. Chem. A 2009, 113, 2075. (41) Murdachaew, G.; Mundy, C. J.; Schenter, G. K. J. Chem. Phys. 2010, 132, 164102. (42) Guillot, B. J. Mol. Liq. 2002, 101, 219. (43) Ball, P. Chem. Rev. 2008, 108, 74. (44) See several papers in Water: A Comprehensive Treatise; Franks, F., Ed.; Plenum: New York, 1972; Vol. 1. (45) The CP2K Developers Group. http://cp2k.berlios.de/ (20002010). (46) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford: New York, 1987. (47) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, CA, 1996. (48) Aguado, A.; Madden, P. A. J. Chem. Phys. 2003, 119, 7471. (49) Laino, T.; Hutter, J. J. Chem. Phys. 2008, 129, 074102. (50) Klopman, G. J. Am. Chem. Soc. 1964, 86, 4550. (51) Dewar, M. J. S.; Hojvat, N. L. J. Chem. Phys. 1961, 34, 1232. (52) Dewar, M. J. S.; Hojvat, N. L. Proc. R. Soc. London, Ser. A 1961, 264, 431. (53) Dewar, M. J. S.; Sabelli, N. L. J. Phys. Chem. 1962, 66, 2310. (54) Ohno, K. Theor. Chim. Acta 1964, 2, 219. (55) Kuo, I.-F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I. J. Chem. Theory Comput. 2006, 2, 1274.

ARTICLE

(56) Nose, S. J. Chem. Phys. 1984, 81, 511. (57) Nose, S. Mol. Phys. 1984, 52, 255. (58) Hoover, W. G. Phys. Rev. A 1986, 34, 2499. (59) Nose, S. Mol. Phys. 1986, 57, 187. (60) Martyna, G. J.; Klein, M. L.; Tuckerman, M. E. J. Chem. Phys. 1992, 97, 2635. (61) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. J. Chem. Phys. 1993, 99, 2796. (62) Kuo, I.-F. W.; Mundy, C. J. Science 2004, 303, 658. (63) Kuo, I.-F. W.; Mundy, C. J.; Eggimann, B. L.; McGrath, M. J.; Siepmann, J. I.; Chen, B.; Vieceli, J.; Tobias, D. J. J. Phys. Chem. B 2006, 110, 3738. (64) Clough, S. A.; Beers, Y.; Klein, G. P.; Rothman, L. S. J. Chem. Phys. 1973, 59, 2254. (65) Tu, Y.; Laaksonen, A. Chem. Phys. Lett. 2000, 329, 283. (66) Silvestrelli, P. L.; Parrinello, M. Phys. Rev. Lett. 1999, 82, 3308. (67) Badyal, Y. S.; Saboungi, M.-L.; Price, D. L.; Shastri, S. D.; Haeffner, D. R.; Soper, A. K. J. Chem. Phys. 2000, 112, 9206. (68) (a) Soper, A. K. Chem. Phys. 2000, 258, 121. See also: (b) Disordered Materials Website. http://www.isis.stfc.ac.uk/groups/ disordered-materials/disordered-materials-6252.html (2011). (69) Paesani, F.; Iuchi, S.; Voth, G. A. J. Chem. Phys. 2007, 127, 074506. (70) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. J. Chem. Phys. 2008, 128, 094314. (71) Mills, R. J. Phys. Chem. 1973, 77, 685. (72) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (73) Yeh, I.-C.; Hummer, G. J. Phys. Chem. B 2004, 108, 15873. (74) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2001, 114, 363. (75) Dang, L. X.; Chang, T.-M. J. Chem. Phys. 1997, 106, 8149. (76) Soper, A. K.; Bruni, F.; Ricci, M. A. J. Chem. Phys. 1997, 106, 247. (77) Kumar, R.; Schmidt, J. R.; Skinner, J. L. J. Chem. Phys. 2007, 126, 204107. (78) Properties of Ordinary Water-Substance; Dorsey, N. E., Ed.; Hafner: New York, 1968. (79) Hura, G.; Sorenson, J. M.; Glaeser, R. M.; Head-Gordon, T. J. Chem. Phys. 2000, 113, 9140. (80) Sorenson, J. M.; Hura, G.; Glaeser, R. M.; Head-Gordon, T. J. Chem. Phys. 2000, 113, 9149.

6053

dx.doi.org/10.1021/jp110481m |J. Phys. Chem. A 2011, 115, 6046–6053