QMPFF2 Performance in Gas Phase. QMPFF2 parameters were fitted to QM data on properties of 144 molecules and 79 molecular dimers, the total number of dimer conformations being 2,916. The relative rms deviations (RMSD) for absolute values of dipole moments and quadrupole and polarizability tensors were 6.2%, 23% and 3.7%, respectively. QMPFF2 fits the QM energies of all of the dimers in the training set well with a RMSD of 0.38 kcal/mol. For the 165 training dimers that have been optimized to have minimum energy, the RMSD of the energy and geometry calculated by QMPFF2 relative to ab initio QM were 0.50 kcal/mol and 0.09 Å.
Table 1 compares QMPFF2 and QM values for basic properties of the water molecule: absolute value of dipole moment, eigenvalues of quadrupole, and polarizability tensors. As seen, QMPFF2 simulates dipole moment exactly, whereas relative errors are 5.0% and 3.2% for quadrupole and polarizability, respectively.
As for interactions of water molecules, Fig. 6, which is published as supporting information on the PNAS web site, illustrates the typical accuracy of QMPFF2 fit for different conformation sets of water dimer, whereas Table 2 compares QMPFF2 and QM energies for optimized water multimers. As seen from Table 2, for the optimized dimer conformations (presented in the QMPFF2 training set) the overall accuracy is 0.4 kcal/mol, which is comparable with the accuracy for dimers of other molecules. On the other hand, QMPFF2 performance outside the training set is validated by comparison of calculated and QM energies of water multimers, which is a challenging test because of the role of many-body effects. As seen the agreement (including even correct ranking of the hexamer energies) is excellent considering that no systems larger than dimer were used in the parameterization of QMPFF2.
QMPFF2 also reasonably predicts dimerization Gibbs energies: rms deviations between calculated and experimental values was found to be 0.35 kcal/mol for 27 different homodimers and heterodimers, consistent with the overall QMPFF2 accuracy. These values, which are closely related to the second virial coefficients, were calculated by using standard methodology neglecting the molecule flexibility with the quantum effects being taken into account in the framework of semiclassical approach. However, we have found from careful path integral calculations that the molecular flexibility may contribute rather essentially to the virial coefficient values (up to several dozen percent in the relevant temperature range), especially in the case of water vapor for which the quantum effects are not negligible.
Bulk Properties. Classical and PIMD simulations have been performed to calculate water characteristics in liquid phase (see Materials and Methods for details of MD protocols). The bulk water properties calculated from the MD trajectories are compared in Table 3 with experimental data taken from ref. 14, except for the diffusion constant taken from ref. 15. As seen, the QMPFF2 classical MD results are close to the experimental values. With quantum MD the agreement is even better. Quantum MD decreases the liquid binding energy by 1.3 kcal/mol to give a value within 3% of experimental. A similar improvement is seen for the heat capacity. The temperature dependence of the liquid binding energy is also reasonably described by QMPFF2 as seen in Fig. 1.
Fig. 1 indicates that the classical simulations resulted in essential underestimation of the liquid binding energy, whereas the quantum effects taken into account by PIMD calculations contribute up to 1–1.5 kcal/mol within the considered temperature range, bringing the curve much closer to the experimental one. This contribution is caused mainly by intermolecular motions, which are adequately treated by calculations with a comparatively low PIMD discretization index, P = 4. As for intramolecular vibrations, the change of QM zero point energy transferring from gas to liquid is small for harmonic intramolecular potentials used in QMPFF2 so the corresponding corrections can be neglected. On the other hand, for essentially anharmonic potentials the correction could lower the liquid binding energy by about a half of kcal/mol (14). Note also that the classical simulations with rigid molecules result in less-bonded liquid phase, which mimics the quantum effects so the liquid binding energy is incidentally very close to the PIMD result (see Fig. 1).
It follows from Table 3 that the greatest improvement provided by quantum MD is a >50% increase in the diffusion constant [which is in agreement with other results (18 and 19)] to fit experiment much better. This finding is also illustrated by Fig. 2 where the temperature dependence of water diffusivity is plotted as found from both classical and PIMD calculations. Note that this QM correction can only be considered an approximation, as the traditional path integral technique does not simulate kinetic properties with sufficient accuracy.
Water Has a Density Maximum. At room temperature, the QMPFF2 density calculated with either classical or quantum MD fits experiment to 0.5%. The QMPFF2 results for other temperatures are plotted in Fig. 3. As seen, the classical MD curve reproduces the density maximum, which is certainly the best-known anomalous property of water. The calculated value ρmax = 1.005 g·cm−3 is in very good agreement with the experimental value 1.000 g·cm−3. The temperature Tmax of the maximum density turns out to be near 8°C, which is between the experimental result 4°C and the value of ≈14°C recommended in ref. 20 for comparison with the results of classical MD. More generally, QMPFF2 with classical MD reproduces the density of water over the entire range −25°C to 100°C to within 0.6% (Fig. 3). With quantum MD, the results are even slightly better for temperatures >0°C, but the density for lower temperatures is incorrectly high. The reason for this discrepancy probably results from problems with the MD algorithm in the supercooled region but more investigation is needed.
Radial Distribution Functions. The radial distribution functions that characterize the structure of liquid water are presented in Fig. 4 for both classical and quantum MD. The QMPFF2 positions of all of the extremes agree very well with the experimental values to within 0.1 Å. The amplitudes are fairly good for OO and HH, but the first OH peak height is overestimated, which is also typical for other models. The quantum MD simulations result in less structured water, giving better overall description of the experiment than classical MD, in agreement with other results (14, 18, 19). Thus, the amplitude of the first OO peak is decreased from the classical MD value of 3.17 to the quantum MD estimate of 2.87, which approaches the experimental value (21, 22) of 2.8. Also, the height of the first OH peak is decreased from 1.58 to 1.30, which should be compared with the experimental value (21) of 1.10. As for the first HH peak, quantum result is worse than the classical, which is likely because of too soft flap angle bending in potential energy surface for water dimer.
Comparison to Specialized Water Models. The quality of QMPFF2 predictions of water properties is overall at least as good as the best specially designed water models. For example, the classical and QM estimates of the liquid binding energy are respectively −11.34 and −9.8 kcal/mol for the MCDHO potential (14); for other ab initio models the classical energy is −10.65, −9.49, and −11.2 kcal/mol for NCC (10), NEMO (11), and TTM2-R (13), respectively. The diffusion coefficient at ambient conditions is found to be 2.6·10−5, 1.3·10−5, and 2.2·10−5 cm2·s−1 for the NCC, NEMO, and TTM2-R potentials. As these values have been determined from classical MD simulations, the diffusion constants should probably be increased by a factor ≈1.5 or greater to account for the increased mobility with quantum MD. Hence, the values for TTM2-R and NCC are less accurate than appears at first glance and the value for NEMO is more accurate.
Most significant is the success of QMPFF2 in reproducing the density/temperature relationship. Water density at room temperature was found to be 1.02 ± 0.01, 0.983, and 1.046 g·cm−3 for MCDHO, NEMO, and TTM2-R, respectively (11–13), which differ from the experimental value by 2–4%, compared with 0.5% for QMPFF2. Moreover, for these models, there are no data on the temperature dependence of the density except for NEMO, where no density maximum was found. On the other hand, a number of empirical water models do find a density maximum, with the Tmax value varying from −40°C to +30°C; nevertheless, most of these models give too sharp a dependence of density on temperature (20) (e.g., see the curve for one of the best models, TIP5P, plotted in Fig. 3).