We report here a new force field for water based solely on quantum mechanics (QM) calculations with no empirical data. The QM was at a high level, coupled cluster single double triple, for all orientations and distances for water dimer plus X3LYP density functional theory (DFT) on 19 larger water clusters. In addition, we included charge and polarization based on the polarizable charge equilibration method and nonbond interactions from DFT-D3 calculations on the H2 and O2 crystal. This model, denoted as RexPoN, provides quite excellent agreement with experimental (expr) data for the solid and liquid phase of water: Tmelt= 273.3 K (expr = 273.15 K) and properties at 298 K: ΔHvap = 10.36 kcal/mol (expr = 10.52), density = 0.9965 gr/cm3 (expr = 0.9965), entropy = 68.4 (J/mol)/K (expr = 69.9), dielectric constant = 76.1 (expr = 78.4), and ln Ds (self-diffusion coef) = −10.08 (expr = −11.24). Such an accurate force field for water will, we believe, be useful for full solvent calculations of electrocatalysis, where we can restrict QM water to just the first one or two layers involving reactions, using RexPoN to provide the polarization for a more distant solvent. Also, RexPoN may provide a better description of the solvent for proteins, DNA, polymers, and inorganic systems for applications to biomolecular, pharma, electrocatalysis (fuel cells and water splitting), and batteries where interaction with explicit water molecules plays a significant role.

Because of its central role in life and sustenance, the properties of water have been widely studied. Although chemically simple, it has been hard to model1,2 with many anomalies such as the critical point of supercooled water below 219 K,3 the maximum density at 4 °C, the high dielectric constant and boiling point, the non-monotonic isothermal compressibility and specific heat,4 and the fifteen known crystal structures for ice.5,6 Because of the lack of long range order, molecular dynamics (MD) simulations have played the major role in elucidating these properties. These MD studies are nearly all based on empirical force fields (FFs) because quantum mechanics (QM) simulations are limited to 100s of atoms for 10s of picoseconds (ps). Starting with the first model of Bernal and Fowler in 1933 and the first computer simulation of water by Barker and Watts in 1969, hundreds of water models have been developed and applied to the properties of water1,2 but there remain many puzzles.

The most popular empirical FFs1,2 are the TIPnP (n = 3, 4, 5)7–9 and simple point-charge (SPC)10 families, which use rigid molecular structures. The parameters in these models are fitted to reproduce one or more macroscopic properties of water but none reproduce all important properties of the solid and liquid phase.

Flexible (e.g., TIP4PF11 and SPC/Fw12) and fluctuating charge (e.g., TIP4PQ77 and SPC/FQ13) models are used less often because of increased costs with no obvious increase in accuracy.

Polarizable FFs [e.g., Anisotropic Site Potentials (ASP);14 Symmetry-adapted Perturbation Theory (SAPT);15 and Thole-type Model (TTM) families16] and Drude oscillator methods (e.g., SWM4-DP17) have been used a little, but some exhibit a “polarization catastrophe” at short distances.18 

Some water models have been developed from fitting many body potential energy terms to millions of high-quality energy datasets [e.g., Huang, Braams, and Bowman (HBB);19 CC-pol;20 and MB-pol21]. These models are not easily transferable to other systems since for any new set of interactions (new species), this tedious fitting procedure must be repeated.22 

Ab initio molecular dynamics (AIMD) methods have been improved over the last two decades to perform at least close to empirical models.23 Various computational techniques have been implemented to improve the qualities of water models but not yet successful. For example, machine learning methods have been used to improve the quality of explicit polarization and Coulombic terms in the water models.24 More discussion of different types of water models is provided elsewhere.1,2

We decided to develop a FF based fully on the best available QM, with the hope that it would correctly describe the standard properties of water including the melting point (Tm), density (ρ), heat of vaporization (ΔHvap), entropy (S0), dielectric constant (ε), self-diffusion constant (Ds), oxygen-oxygen radial distribution function (ROO), and which might additionally explain the anomalies of water. Our procedure is as follows:

  • First, establish the long range interactions: the van der Waals nonbond (NB) interactions (London dispersion plus Pauli repulsion) from QM.

  • Describe the electrostatic (polarization and charge) interactions based on our recently developed polarizable charge equilibration (PQEq) method that has been shown to reproduce the QM polarization energies due to bringing point dipole up to various small molecules.

  • Use Coupled Cluster Single Double Triple, CCSD(T), QM on a water monomer25 (including OH bond dissociation and HOH angle terms) to determine the bond and angle valence terms.

  • Use Coupled Cluster Single Double Triple, CCSD(T), QM on a water dimer25 to determine the hydrogen bond (HB) corrections.

The final FF is called RexPoN to indicate that it includes charge and polarization (Po) and Nonbond (N) while allowing bond breaking (Rex). RexPoN is aimed at doing both reactive and nonreactive systems. Thus, we use the full RexPoN formalism to define the valence bond and angle terms at the same level as would be used in reactions.

However, in the validations reported here, we use a rigid water model where OH bonds and HOH angles are kept fixed at the values from CCSDT throughout the simulations (see below). There are questions concerning whether the rigid water model describes nuclear quantum effects (NQE) including zero-point energy (ZPE) and tunneling. The hydrogen bonding between water molecules can be affected by ZPE effects for OH bond stretch and HOH bending vibrational modes.26 However, it has been shown that previous estimates of NQE effects were greatly overestimated because of the use of harmonic intramolecular interactions, ignoring anharmonicity.27 In addition, proton delocalization and tunneling effects are expected to be important only at very low temperatures and high pressures which are not the conditions of simulations here.26 Moreover, RexPoN includes the correct anharmonicity at the level of CCSDT since the OH bond and HOH angle curves were fitted to this level. Therefore, flexible RexPoN will likely do well on NQE effects in water.

The rest of this paper is organized as below. Section II describes the model and the physical motivation behind each potential energy term. Section II also discusses the training sets and optimization process of the FF. Section III provides the results of the simulations for several ice Ih and water properties and compares the results with experimental data. The conclusions are given in Sec. IV.

The total energy of the system is expressed as

Etotal=Enonbond+Evalence,
(1)

where Enonbond is the nonbond (NB) energy and Evalence is the bonded energy terms. Enonbond and Evalence are defined by

Enonbond=Eelect+EvdW+EHB,
(2)
Evalence=Ebond+Eangle,
(3)

where Eelect is the electrostatic, EvdW is the van der Waals (vdW), and EHB is the hydrogen bond (HB) energy. Evalence includes the covalent OH bond (Ebond) and the HOH angle (Eangle) energy terms.

We describe the electrostatic energy of the system using our polarizable charge equilibration (PQEq) model.28,29 In addition to atomic charges, we use PQEq to include the polarization arising from interatomic interactions and external electric fields. Each atom is considered to have a core including all the mass of the atom and an electron shell that is massless. The core and shell are described by 1s Gaussian shaped electron densities (ρ) with a finite distribution of charge on each atom rather than fixed point charges used by many other FFs (see Fig. 1). The core has a variable charge (qi) and a fixed charge (+1) part. The shell only has a fixed charge of −1.The dynamic position of the shell is found by balancing the external forces due to other atoms and applied fields with the spring force constant (Ks) that connects the shell and core of the atom. The electrostatic interaction energy between the electron densities (ρ) of atoms i and j is given by

Cik,jlr=1rerfαikαjlαik+αjlr,
(4)

where k and l denote the core (c) and shell (s), respectively. The width of the Gaussian distribution (αik and αik) is equal to 0.2341/R2, where R is the core (Rc) or shell (Rs) radius. The electrostatic energy (Eelect) is defined by

Eelectric,ris,qi=iNχi0qi+12Jii0qi2+12Ksric,is2+ik>jlTrik,jlCik,jlrik,jlqikqjl,
(5)

where χi0 is the Mulliken electronegativity, χi0 = (IP + EA)/2, Jii0 is the idempotential, Jii0 = IP − EA, rik,jl is the interatomic distance, IP = atomic ionization potential, and EA = atomic electron affinity. Here T(r) is the 7th order taper function [see Eq. (29)]. The second sum computes the electrostatic energy between the core and shells of all atoms in the system. The charge on the core (qi) is updated every time step via efficient schemes. The massless shell relaxes instantaneously to align itself along the external force vector with no inertial delay. The complete description of the PQEq model and each term in the energy equation is provided in previous publications.28,29 We validated the accuracy of PQEq by comparing with the polarization in QM.28 We reoptimized the PQEq parameters of water for two reasons.

FIG. 1.

The PQEq model for a water molecule. Each atom has a core and shell. The core has a variable (ρH and ρO) and a fixed charge distribution (ρc) both positioned at the center of the atom. The shell has a fixed shell charge distribution (ρs) with negative charge connected by a harmonic spring to ρc. The ρc and ρs have the overall charges of +1 and −1, respectively. The figure shows the positions of shells for a single isolated water molecule. The position of each shell is found by balancing the external forces due to other atoms with the spring force that connects the shell and core of the atom.

FIG. 1.

The PQEq model for a water molecule. Each atom has a core and shell. The core has a variable (ρH and ρO) and a fixed charge distribution (ρc) both positioned at the center of the atom. The shell has a fixed shell charge distribution (ρs) with negative charge connected by a harmonic spring to ρc. The ρc and ρs have the overall charges of +1 and −1, respectively. The figure shows the positions of shells for a single isolated water molecule. The position of each shell is found by balancing the external forces due to other atoms with the spring force that connects the shell and core of the atom.

Close modal

First, to ensure that the charge distribution on oxygen and hydrogen atoms is in reasonable ranges when two water molecules get very close to each other. This is important for describing the equation of state (EOS) of water at high pressures. Thus, we computed the QM [electrostatic potential (ESP) and Mulliken] charge distributions on the water dimer over a range of short distances and adjusted the PQEq parameters such that they produce reasonable charges (see Fig. S1 of the supplementary material).

Second, during the optimization of HB attractive and repulsive terms (see below), we optimize the PQEq parameters such that the electrostatics energy provides the most accurate description of the water dimer energy at long distances. The PQEq parameters of H and O are provided in Table I.

TABLE I.

PQEq parameters for hydrogen and oxygen atoms in a water molecule.

Atomχi0 (eV)Jii0 (eV)Rc= Rs (Å)Ks [(kcal/mol)/Å2]
4.5280 17.9841 0.4452 2037.2006 
8.7410 13.3640 0.8028 814.0445 
Atomχi0 (eV)Jii0 (eV)Rc= Rs (Å)Ks [(kcal/mol)/Å2]
4.5280 17.9841 0.4452 2037.2006 
8.7410 13.3640 0.8028 814.0445 

To determine the vdW nonbond interactions, we calculated the equation of states (EOS) of solid hydrogen and solid oxygen based on the density functional theory (DFT).30,31 For solid hydrogen, we use the Perdew-Burke-Ernzerhof (PBE)32 flavor of DFT, while for a solid oxygen crystal, we use Becke three-parameter Lee-Yang-Parr (B3LYP)33,34 including the hybrid terms important for describing the O2 triplet state35,36 (see below).

1. Solid hydrogen

X-ray diffraction shows that the H2 crystal has a hexagonal closed packed structure with the P63/mmc space group and a c/a ratio close to the ideal value of 1.633.37 However, since hydrogen has no core electrons, X-ray experiments cannot determine the exact position of the H atom coordinates. Therefore, extensive theoretical studies have been carried out to determine the H2 molecule orientation. It was found that the Pca21 space group (see Fig. S2 of the supplementary material) containing a hexagonal close packed structure with four molecules per unit cell is the most stable structure at high pressure ranges (∼110–150 GPa).38,39 Therefore, we selected a Pca21 space group and used the PBE-D3 to compute energy versus volume with single point energy calculations at various volumes. This leads to the PBE-D3 energy curve shown in Fig. 2(a) for densities starting with the volume of 17.440 Å3 (ρ = 2.626 cm3/mol) corresponding to the pressure of 120 GPa at 300 K and passing through the volume of 44.082 Å3 (6.627 cm3/mol) corresponding to the pressure of 10 GPa at 300 K. We also scale the simulation cell to large volumes as high as 424.572 Å3 (63.920 cm3/mol) to capture the long-range interactions. The HH bond distance was fixed at the equilibrium value of 0.74 Å for all volumes to avoid any contribution to the total PBE-D3 curve from changes in the HH bond length. All calculations in this section were performed using the Vienna Ab initio simulation package (VASP) (version 5.3.5).40 The generalized gradient approximation (GGA) was used for the exchange–correlation energy32 with a plane wave energy cutoff of 600 eV. We used the convergence criteria of 10−6 eV for energy and 10−4 eV Å−1 for force. The k-space sampling was the gamma point centered with full space group symmetry.

FIG. 2.

The change of vdW energies with volume for (a) solid hydrogen and (b) solid oxygen computed by DFT-D3 (open circles) and RexPoN (red squares). The vdW curves are the sum of the nonbond repulsion (NBrep) and nonbond attraction (NBatt). The dotted vertical line in the insets shows the volume for the minimum energy. The DFT-D3 energy in (a) is computed by PBE-D3, and the volumes of 17.440 and 44.082 Å3 are related to the experimental pressures of 120 and 10 GPa at 300 K, respectively. The DFT-D3 energy in (b) is computed by B3LYP-D3, and the volumes of 17.369 and 41.376 Å3 are related to the experimental pressures of 300 and 10 GPa at 300 K, respectively.

FIG. 2.

The change of vdW energies with volume for (a) solid hydrogen and (b) solid oxygen computed by DFT-D3 (open circles) and RexPoN (red squares). The vdW curves are the sum of the nonbond repulsion (NBrep) and nonbond attraction (NBatt). The dotted vertical line in the insets shows the volume for the minimum energy. The DFT-D3 energy in (a) is computed by PBE-D3, and the volumes of 17.440 and 44.082 Å3 are related to the experimental pressures of 120 and 10 GPa at 300 K, respectively. The DFT-D3 energy in (b) is computed by B3LYP-D3, and the volumes of 17.369 and 41.376 Å3 are related to the experimental pressures of 300 and 10 GPa at 300 K, respectively.

Close modal

2. Solid oxygen

For solid oxygen, we selected the stable α-phase crystal structure, which according to neutron-diffraction and X-ray measurements has a monoclinic base-center structure with a C2/m space group35,36 (see Fig. S2 of the supplementary material). For DFT, we used the B3LYP hybrid method. This is because O2 has a 3Σg triplet ground state that is poorly described by PBE. The ground state of O2 crystal has the two O2 molecules coupled antiferromagnetically to obtain a singlet spin state. Since the DFT formalism uses a single Slater determinant, it cannot properly describe the single spin state, often leading to severe convergence problems. Instead, we describe the O2 crystal assuming a high spin state quintet state (S = 2, with 4 extra up spin orbitals per cell, MS = 2) which is correctly described with a single Slater determinant using a spin unrestricted DFT formalism. Our computed heat of vaporization of solid oxygen using B3LYP-D3 and high spin O2 is −2.3 kcal/mol after correcting for zero-point energy, which is very close to the experimental value of −2.1 kcal/mol.36 The EOS of solid oxygen was then obtained by performing single point energy calculation at different volumes. This leads to the B3LYP-D3 energy curve shown in Fig. 2(b). The smallest volume of 17.369 Å3 (ρ = 5.230 cm3/mol) corresponds to a pressure of 300 GPa at 300 K,41,42 the volume of 41.376 Å3 (12.459 cm3/mol) corresponds to a pressure of 10 GPa at 10 K,43 and the volume 69.439 Å3 (20.909 cm3/mol) corresponds to a pressure of 0 GPa at 10 K.43 The OO bond distance was fixed at the equilibrium value of 1.207 Å. All calculations in this section were performed using the CRYSTAL package (version 14)44 using the modified Gaussian basis set, m-6-311G(2d) for oxygen.45 We used the convergence criteria of 10−7 a.u. for energy. The k-space sampling was the gamma point centered with full space group symmetry.

We used the CRYSTAL package for solid oxygen because the VASP package is over 1000 times slower than CRYSTAL for B3LYP calculations. However, we used VASP for PBE-D3 calculations of solid hydrogen because it is far easier to use than CRYSTAL and is only a factor of 3 slower.

The hydrogen and oxygen DFT-D3 energy curves were then fitted to obtain the first principles based vdW potential energy curves for H and O. This vdW energy is the sum of a monotonic attractive nonbond term (NBatt) and a monotonic repulsive non-bond term (NBrep). NBatt has the form of −C6/r6 which means that at sufficiently long r the wavefunction on separate atoms does not overlap. The NBatt must remain finite as the distances become sufficiently small. Therefore, we define NBatt energy (ENBatt) as

ENBattr=i<jC6ijrij6+RvdWij6,
(6)

which is similar to low gradient formulation46,47 that was also used in ReaxFF-lg.48 Here, r is the interatomic distance, C6 is the dispersion energy correction parameter, and RvdW is the equilibrium vdW distance between atoms i and j. The NBrep energy (ENBrep) is purely repulsive (positive) and monotonic which we expand using a group of exponential functions as

ENBrepr=i<jAijexpαijrij+βijexp(γijrijnij+ηijrij+δij),
(7)

where A, α, β, γ, n, η, and δ are the parameters.

In developing RexPoN, we want to use generic approaches likely to have the correct systematics as a function of the rows and columns of the periodic table. Thus we used as initial values for RvdW and C6 the values from the universal force field,49 which has parameters up to Lw (Z = 103). We expect that the RvdW should be reasonably accurate; however, the C6 was based on calculated atomic polarizabilities and needs to be adjusted. Next, we define the initial parameters of ENBrep knowing that ENBrep(r) = EvdW(r)-ENBatt(r). Then, the parameters of ENBatt and ENBrep are optimized to provide th6e best fit to the Evdw curve.

To describe HO vdW interaction, we used the geometric mean average of the two-body vdW potential energy curves of H and O as in the following equations:

ENBreprOH=ENBreprHHENBrep(rOO),
(8)
RvdW,OH=RvdW,HHRvdW,OO,
(9)
C6,OH=C6,HHC6,OO.
(10)

The final ENBatt and ENBrep parameters are given in Tables II and III, respectively.

TABLE II.

ENBatt parameters of the vdW energy term of the water model.

Atom iAtom jRe (Å)C6 [(kcal Å6)/mol]
3.2692 239.4803 
2.9938 634.7066 
3.1285 389.8714 
Atom iAtom jRe (Å)C6 [(kcal Å6)/mol]
3.2692 239.4803 
2.9938 634.7066 
3.1285 389.8714 
TABLE III.

ENBrep parameters of the vdW energy term of the water model.

Atom iAtom jA (kcal/mol)αβγnηδ
0.0212 −2.4002 10.5460 −0.8113 0.0014 −0.2607 −0.5640 
0.0272 −3.6468 14.5278 −0.3235 0.0305 0.0331 −0.5688 
0.0906 −3.1455 10.0823 −0.0007 −1.0000 0.0000 0.0000 
Atom iAtom jA (kcal/mol)αβγnηδ
0.0212 −2.4002 10.5460 −0.8113 0.0014 −0.2607 −0.5640 
0.0272 −3.6468 14.5278 −0.3235 0.0305 0.0331 −0.5688 
0.0906 −3.1455 10.0823 −0.0007 −1.0000 0.0000 0.0000 

3. Validation of the method

We tested the accuracy of DFT-D3 by comparing to the experimental pressure-volume (P-V) EOS of solid hydrogen. We used experimental measurements on solid hydrogen isotopes (H2 and D2) within the pressure range of 10–120 GPa at 300 K for these comparisons37 (Fig. 3). We first minimized the atomic positions related to this pressure range (i.e., from 17.440 to 44.082 Å3) using PBE-D3 to obtain the total energy at 0 K. Then, we corrected this energy for the zero point energy (ZPE) contributions and also free energy (thermal effects) at 300 K (G = H-TS). The total PBE-D3 energy before and after correction at 300 K is shown in Fig. S4 of the supplementary material. Then, the total energy curve at 300 K was fitted to third-order Birch-Murnaghan isothermal EOS,50EBM(V). The P-V EOS was then obtained by

PBMV=EBMVS.
(11)

As shown in Fig. 3, our PBE-D3 predicted EOS is in good agreement with experiment, indicating that DFT-D3 is sufficiently accurate to describe the vdW interactions of dense solid systems.

FIG. 3.

The computed and experimental pressure-volume equation of state of the hydrogen. The DFT curve (black) has been corrected for zero-point energy and thermal effects at 300 K (red) to compare with the experimental curve37 (blue) at 300 K.

FIG. 3.

The computed and experimental pressure-volume equation of state of the hydrogen. The DFT curve (black) has been corrected for zero-point energy and thermal effects at 300 K (red) to compare with the experimental curve37 (blue) at 300 K.

Close modal

1. The bond order (BO) expression

To describe the covalent O–H bond, we used the bond order (BO) concept to capture the dependence of the bonding on the overlap (actually exchange kinetic energy) between orbitals on the bonding atoms. Here we use a sigmoid function as in ReaxFF51 to define the bond distance to the bond order relation

BOij=1.101.0+exppbo1rijpbo2Repbo2,
(12)

where pbo1 and pbo2 are the parameters, Re is the equilibrium bond distance, and r is the bond distance between atom i and atom j. We choose the BO = 1.0 at Re which results in pbo1 = 2.3026. This allows the formal bond order to approach 1.1 for very small values of r. Therefore, the only free parameter that needs to be determined is pbo2. To make an initial good guess for pbo2, we choose the BO = 0.5 at the inflection point distance (RH/2) of the bond energy curve. The QM bond energy curves (see below) are used to obtain pbo2 parameters for H2, O2, and OH (in H2O molecule) cases. The final bond order curves are shown in Fig. S5 of the supplementary material, and the BO parameters are summarized in Table IV.

TABLE IV.

Bond order parameters for HH, OO, and OH bonds.

Bond typeRe(Å)pbo1pbo2
H–H 0.7450 2.3026 1.3532 
O–O 1.2070 2.3026 1.6150 
O–H 0.9587 2.3026 1.4428 
Bond typeRe(Å)pbo1pbo2
H–H 0.7450 2.3026 1.3532 
O–O 1.2070 2.3026 1.6150 
O–H 0.9587 2.3026 1.4428 

We define the bond energy (Ebond) to match the bond dissociation curves in H2, O2, and OH (in H2O) molecules for the highest quality QM calculations. Thus, given Evdw and Eelect, we define EQ2 and Ebond as

EQ2rij=EQMriji,jNNEvdWrkl+NEelectrkl,
(13)
Ebondrij=EQ2rijEvdWrij,
(14)

where EQM is the total QM energy due to the stretching bond between atom i and j. k and l are the atom indices (k, li, j). EQ2 is the net two-body QM energy between atom i and j, which is obtained by subtracting from EQM the EvdW between all atom pairs except for i and j and also subtracting the total Eelect of the system [Eq. (5)]. Therefore, based on EvdW and Eelectro which are already established, the Ebond is defined automatically to match the QM bonding energy curves.

We used the exact ab initio singlet energy data for the HH bond scan in a H2 molecule52 and high quality ab initio and spectroscopic energy data for the OO bond scan in an O2 molecule.53,54 For diatomic H2 and O2 molecules, the Evdw(rkl) and Eelecto(rkl) sums in Eq. (13) are zero and therefore EQ2(rij) = EQM(rij). For the OH bond in the H2O molecule, we used the accurate water monomer potential energy surfaces (Bowman PES) obtained from CCSD(T) calculations.25 Here, we started from the equilibrium water structure (OH = 0.9572 Å and HÔH = 104.52°) and stretched one of the OH bonds to infinity while the angle and other bond are fixed. Then, the EQM was computed at each distance using Bowman PES. Here, the HH Evdw and total Eelectro [Eq. (13)] should be subtracted from EQM to get the net OH two-body curve (EQ2). The EQ2 curve is fitted to an extended-Rydberg type function of the form

ERydrij=D0expb0a0rRec0+b0a0rRe×c0+a0rRec0+a0(rRe),
(15)

where D0, a0, b0, and c0 are the parameters to be determined. Figure 4 shows the energy components and the fit to the EQ2 curve for H2, O2, and OH bond dissociation. The final parameters of Eq. (15) are given in Table V.

FIG. 4.

The QM (black circles) and RexPoN (red) total energies for bond dissociations of (a) HH, (b) OO, and (c) OH. The energy components of RexPoN are also shown. Insets in (b) and (c) show the curves for a broader range such that the complete range of the bond energy curve (blue) is shown.

FIG. 4.

The QM (black circles) and RexPoN (red) total energies for bond dissociations of (a) HH, (b) OO, and (c) OH. The energy components of RexPoN are also shown. Insets in (b) and (c) show the curves for a broader range such that the complete range of the bond energy curve (blue) is shown.

Close modal
TABLE V.

The parameters of the ERyd(r) bonding energy term.

Bond typeD0 (kcal/mol)Re (Å)a0b0c0
H–H 109.4853 0.7450 1.0000 4.0000 1.0000 
O–O 120.4881 1.2070 1.0000 5.2000 1.0000 
O–H 113.2500 0.9587 0.5577 7.1812 1.0000 
Bond typeD0 (kcal/mol)Re (Å)a0b0c0
H–H 109.4853 0.7450 1.0000 4.0000 1.0000 
O–O 120.4881 1.2070 1.0000 5.2000 1.0000 
O–H 113.2500 0.9587 0.5577 7.1812 1.0000 

To define the angle term, we fit a simple cosine angle function to the Bowman PES for the water monomer,25 

Eanglerij,rjk,θijk=2Csin2(θ0)F1BOijF2BOjk×cos(θijk)cos(θ0)2,
(16)
F1BOij=1expaBOijb,
(17)
F2BOjk=1expaBOjkb,
(18)

where C, θ0, a, and b are the adjustable parameters. θijk is the angle between atoms i, j, and k. F1 and F2 are the scaling functions. The scaling functions go to zero as OH gets very large and become equal to one at the equilibrium bond length (rij= Re or rjk= Re). The QM energy as a function of θHOH with both bonds fixed at Re = 0.9572 Å (i.e., F1= F2 = 1) is shown in Fig. 5. The angle term together with Evdw and Eelectro is fitted to the QM curve. The parameters of angle energy term are given in Table VI.

FIG. 5.

The energy components and the fit to the HOH angle scan in the water molecule.

FIG. 5.

The energy components and the fit to the HOH angle scan in the water molecule.

Close modal
TABLE VI.

The parameters of the angle energy (Eangle) term.

Angle typeC (kcal/mol)θ0 (degrees)ab
H–O–H 31.6132 98.8511 3.3562 3.8453 
Angle typeC (kcal/mol)θ0 (degrees)ab
H–O–H 31.6132 98.8511 3.3562 3.8453 

The water dimer presents a deep minimum (∼5.0 kcal/mol) at its global minimum which leads to ice structures having tetrahedral arrangements leading to four strong HBs to each water. The HB involves attraction between a partially positive hydrogen of the donor water molecule to the partially negative lone pair of the acceptor water molecule. Similarly liquid water forms an extensive HB network that is responsible for many of its unique properties. This HB depends strongly on the alignment of the donor hydrogen with respect to the lone pair of the acceptor molecule. In addition, the lone pairs on each water molecule repel the lone-pair electrons of nearby water molecules. To describe their roles, we include the role of lone pairs on each water molecule explicitly in our FF.

1. Explicit lone pairs

Although, the lone pair electrons do not present distinct directed electron densities, there are minima in the electrostatic potential surfaces around the isolated water molecule that suggest a tetrahedral structure for the water molecule.55 The accurate Bowman PES for a water dimer also confirms two local minima (∼109.5° apart) during the rotation of donor water around the acceptor water molecule. A tetrahedral bonding network is also found in condensed phases of water such as hexagonal ice. Therefore, we define tetrahedrally arranged lone pairs (sp3 hybridized with an angle of 109.5° between them) defined with respect to the two OH bonds.

The explicit lone pairs of the acceptor molecule in a water dimer are shown in Fig. 6. The lone-pair unit vectors (uα and uβ) are based on the OH bond unit vectors (u1 and u2) and the θL = 109.5° angle between them. First, the u1 and u2 vectors are used to define the perpendicular (up) and bisector (ub) unit vectors. The up vector is perpendicular to the plane of H2O, up=(u1×u2)/u1×u2, and ub is in the plane of H2O and bisects the HOH angle, ub=u1+u2/2.Then, the vector that connects acceptor oxygen to the hydrogen of the donor is normalized to give the uOH unit vector. The α and β angles (shown in Fig. 6) are defined as the angles that uOH constructs with uα and uβ unit vectors, respectively. Finally, we use the following three conditions to completely define lone-pair vectors (uα and uβ):

  • the angle between lone-pair vectors and ub is θL/2 = 54.75° because ub bisects the angle between uα and uβ,

  • the angle between lone-pair vectors and up is π/2−θL/2 = 35.25°, and

  • lone-pair, ub, and up vectors are all in the same plane.

Equations (19)–(21) are written based on these conditions. They are solved simultaneously to obtain uα and uβ,

ubuα=cosθL2,
(19)
upuα=cosπθL2,
(20)
uαub×up=0.
(21)

The EHB is defined as the sum of attractive (EHBatt) and repulsive (EHBrep) terms

EHBr,θ,α,β=EHBattrkj,α,β+EHBreprki,θijk.
(22)

The EHBatt term is defined based on lone-pair angles (α and β) and the distance (rkj) between the acceptor oxygen (k) and the hydrogen of the donor molecule (j). The EHBatt is expressed as

EHBattrkj,α,β=EHBattrkjEHBattα,β,
(23)
EHBattrkj=hb1hb2e2hb3rkjhb4hb52,
(24)
EHBattα,β=hb61+ehb7sin6α2+sin6β2+hb8,
(25)

where hb1 to hb8 are the adjustable parameters. The EHBrep term is defined by

EHBreprik,θijk=EHBreprikEHBrepθijk,
(26)
EHBreprik=lp1explp2riklp3,
(27)
EHBrep(θijk)=1+lp4sin2θijk2.
(28)

Here, lp1 to lp6 are the adjustable parameters and rik is the oxygen-oxygen distance. We define angle θijk to be between atoms i, j, and k atoms (see Fig. 6) as the average position of the lone pairs is positioned on the center of the oxygen atoms.55 Since the lone pairs are only defined by vectors, we transfer their force components to the hydrogens of the same molecule via chain derivatives.

FIG. 6.

Schematic defining the explicit lone pairs for the acceptor molecule in the water dimer. The lone-pair vectors (uα and uβ) are constructed based on the OH bonds vectors (u1 and u2) and by solving Eqs. (19)–(21). Note that the angle between the LPs was selected as 109.5° without optimization.

FIG. 6.

Schematic defining the explicit lone pairs for the acceptor molecule in the water dimer. The lone-pair vectors (uα and uβ) are constructed based on the OH bonds vectors (u1 and u2) and by solving Eqs. (19)–(21). Note that the angle between the LPs was selected as 109.5° without optimization.

Close modal

For all nonbond energy terms (i.e., Eelectro, EvdW, and EHB), we use a 7th order taper function to cutoff the forces and energies smoothly at the cutoff distance (matching energy and the first three derivatives at r = 0 and r = cutoff). It is given by

Trij=α=07Tapαrijrcutα,
(29)

where rcut is a cutoff distance and Tap7 = 20, Tap6 = −70, Tap5 = 84, Tap4 = −35, Tap3 = 0, Tap2 = 0, Tap1 = 0, and Tap0 = 1. We used a cutoff distance of 12.0 Å for EvdW and Eelectro. W used a cutoff of 4.6 Å for EHB to include only the first nearest neighbors in hydrogen bonding.

2. HB training set

To obtain EHB parameters, we started with a minimal training set scanning only four dimer cases. These four cases were selected to establish the range of the parameters and to avoid overfitting. The QM energy for all these cases was obtained using the Bowman PES. The dimer scans are shown in Fig. 7.

  • Scan 1: one of the water molecules (D2h symmetry) is scanned toward the other molecule in the xy plane. This scan is used to optimize the distance dependent parameters of the EHBrep(rik) term [Eq. (27)]. Here, it is assumed that there is no angle contribution from EHBrep [Eq. (28)] and no contribution from EHBatt.

  • Scan 2: the donor water molecule in the most stable water dimer structure is rotated around the z axis with the oxygen of the donor (Oi) as the center of the rotation.

  • Scan 3: the donor molecule in the most stable water dimer structure is scanned toward the acceptor with the OiHjOk angle equal to 180.0° during the scan. Scan 2 and Scan 3 are used to obtain the EHBrepijk) and EHBatt(rkj) parameters.

  • Scan 4: the donor water molecule in the most stable water dimer is rotating in the xy plane around the oxygen of the acceptor molecule (Ok) as the center of rotation. We used this to optimize the EHBatt(α,β) parameters.

The optimization of the HB parameters for the above four scan cases used conjugate gradients to minimize the root-mean square error (RMSE) of the energy values for each geometry. The energy comparison for the obtained parameter set (FF1) and QM is shown in Fig. 8, which are in very good agreement.

FIG. 7.

Four dimer scan cases used in the training set to optimize HB parameters. In Scan 1, the water molecule on the right is scanned toward the water molecule on the left in the xy plane, in Scan 2, the donor molecule is rotating around the z axis with the Oi as the center of rotation, in Scan 3, the donor molecule is scanned toward the acceptor with OiHjOk = 180.0° and in the xy plane, and in Scan 4, the donor water molecule is rotating in the xy plane around Ok as the center of rotation. In all cases, the bonds and angles are kept fixed during the scans.

FIG. 7.

Four dimer scan cases used in the training set to optimize HB parameters. In Scan 1, the water molecule on the right is scanned toward the water molecule on the left in the xy plane, in Scan 2, the donor molecule is rotating around the z axis with the Oi as the center of rotation, in Scan 3, the donor molecule is scanned toward the acceptor with OiHjOk = 180.0° and in the xy plane, and in Scan 4, the donor water molecule is rotating in the xy plane around Ok as the center of rotation. In all cases, the bonds and angles are kept fixed during the scans.

Close modal
FIG. 8.

Energy comparison between RexPoN and QM using the FF1 (red) and FF2 (blue) optimized parameter sets. The FF1 parameter set is obtained using only the 4 scan cases shown in Fig. 7 while for FF2 in addition to 4 scan cases, the 10 low-lying stationary points of water dimers and water clusters, (H2O)n, containing up to n = 19 clusters were included in the training set.

FIG. 8.

Energy comparison between RexPoN and QM using the FF1 (red) and FF2 (blue) optimized parameter sets. The FF1 parameter set is obtained using only the 4 scan cases shown in Fig. 7 while for FF2 in addition to 4 scan cases, the 10 low-lying stationary points of water dimers and water clusters, (H2O)n, containing up to n = 19 clusters were included in the training set.

Close modal

To improve the quality of the force field for different configurations of water dimer interactions, we also included the 10 low-lying stationary point structures (see Fig. S6 of the supplementary material) of the water dimer in our training set. The reference energies and minimized structures are from the Bowman PES. In addition, to account for possible non-additive HB effects that exist in a bulk system, we expanded our training set by including the global minima of water clusters, (H2O)n, containing up to n = 19 waters. The reference energies and global minimized structures were taken from X3LYP56 hybrid DFT calculations which leads to excellent structures for large water clusters. These calculations were performed using the Jaguar57 software. The cluster structures are provided in Fig. S7 of the supplementary material. Our optimization of the HB term results in a RMSE difference of 0.65 kcal/mol/water for the above 10 water dimers and 19 water clusters. The comparison between RexPoN and QM for the 4 scan cases after including all dimers and water clusters in the optimization (FF2) is shown in Fig. 8. The energy components of RexPoN for FF2 are shown in Fig. S8 of the supplementary material. The systematic lower energy of the FF2 with respect to FF1 is due the non-additive HB effects that are not captured completely by considering just water dimers. Because of that at the OO distance of 2.76 Å, the water dimer energy computed by FF2 is equal to 5.7 kcal/mol which is close to the experimental ice HB energy58 of 5.5 kcal/mol (OO = 2.76 Å). The comparison between the computed energies of RexPoN and QM for water dimers and clusters is provided in Table S1 of the supplementary material. The EHBatt and EHBrep optimized parameters [Eqs. (23)–(28)] for FF2 are provided in Tables VII and VIII, respectively.

TABLE VII.

The parameters of the EHBatt energy term [Eqs. (24) and (25)].

Parametershb1hb2hb3hb4hb5hb6hb7hb8
Value −22.5000 12.7434 1.1345 0.7861 3.5071 1.1693 3.6922 −0.8921 
Parametershb1hb2hb3hb4hb5hb6hb7hb8
Value −22.5000 12.7434 1.1345 0.7861 3.5071 1.1693 3.6922 −0.8921 
TABLE VIII.

The parameters of the EHBrep energy term [Eqs. (27) and (28)].

Parameterslp1lp2lp3lp4
Value 0.7650 10.2511 2.0214 15.0000 
Parameterslp1lp2lp3lp4
Value 0.7650 10.2511 2.0214 15.0000 

To validate the RexPoN FF, we computed several properties of water and ice including the radial distribution function (RDF), heat of vaporization (ΔHvap), dielectric constant (ε), and diffusivity of water at 298 K, and melting point of ice Ih. A quick comparison for several properties of liquid water and also the melting point of ice Ih between experiment, some of the widely used water models, and RexPoN FF is provided in Table IX.

TABLE IX.

Summary of predicted properties from the RexPoN FF compared to experimental data and other popular models.63–65 The melting temperature Tm (K) at 1 atm pressure, standard molar entropy S0 [(J/mol)/K], density ρ (g/cm3), static dielectric constant ε, heat of vaporization ∆Hv (kcal/mol), and self-diffusion coefficient Ds (cm2/s) all at T = 298 K, p = 1 atm. Note that MB-pol is also based only on QM data.

Expt.RexPoNTIP3PTIP4P-2005TIP5PSPC/EMB-pol
Tm 273.15 273.3 146 252 274 215 263.5 
S0 69.9 68.43 72.51 57.47 … 60.30 … 
ρ 0.9965 0.9965 0.98 0.993 0.979 0.994 1.007 
ε 78.4 76.1 94 58 91 68 68.4 
∆Hv 10.52 10.36 10.05 11.99 10.46 11.79 10.93 
ln Ds −11.24 −10.08 −10.2 −11.27 −11.41 −11.08 −10.46 
Expt.RexPoNTIP3PTIP4P-2005TIP5PSPC/EMB-pol
Tm 273.15 273.3 146 252 274 215 263.5 
S0 69.9 68.43 72.51 57.47 … 60.30 … 
ρ 0.9965 0.9965 0.98 0.993 0.979 0.994 1.007 
ε 78.4 76.1 94 58 91 68 68.4 
∆Hv 10.52 10.36 10.05 11.99 10.46 11.79 10.93 
ln Ds −11.24 −10.08 −10.2 −11.27 −11.41 −11.08 −10.46 

We integrated the RexPoN model in the LAMMPS59 MD simulator package and will submit our modifications for inclusion in the official LAMMPS releases. For these water simulations at 298 K, we used a system with 216 molecules per cell, and for ice melting simulations, we used a system with 96 molecules per cell. To control the pressure, we used a barostat with a relaxation time of 1 ps, and to control temperature, we used the Nosé-Hoover thermostat with a damping time of 100 femtosecond (fs).60,61 We used a rigid model for the water molecules since most practical simulations use rigid models. The OH bond distance and HOH angle are kept fixed at their equilibrium values (OH = 0.9572 Å and HÔH = 104.52°) using SHAKE methodology.62 We used a time step of 1.0 fs here, but time steps of 2.0 fs could have been used. We used a time step of 1.0 fs to ensure the accuracy for all properties. We expect that using time steps of 2.0 fs would lead to similar accuracy, but we have not yet tested the precision of computed properties using 2.0 fs time steps.

To determine the melting point of the ice Ih structure, we started with ice Ih at 0 K and increased temperature gradually to the range of 200–320 K during 10 ps followed by 150 ps constant temperature-volume (MD-NVT) calculations. Then, the cell was relaxed for 120 ps during constant temperature-pressure (MD-NPT) simulations. Then at each T, the system was further relaxed for another 100 ps (5 cycles each 20 ps) using MD-NVT for performing the two-phase thermodynamics (2PT) analysis.66,67 The 2PT analysis is used to calculate the thermodynamic properties, including entropy and free energy. 2PT determines the vibrational density of states from the Fourier transform of the velocity autocorrelation function after which corrections are made for the diffusional contributions that lead to a finite density of states at υ = 0. 2PT allows the free energy, entropy, and other thermodynamic properties to be determined from short 20 ps MD trajectories. The accuracy of the 2PT method has been validated against expensive ab initio methods for the calculation of absolute entropy of liquids.64 Figure 9 shows the change of the standard molar entropy (S0) and free energy (A) of the ice versus temperature. Melting is a first-order transition that results in a discontinuity and sudden change in the slope of the S0 and A curves. RexPoN finds that this sharp discontinuity occurs between 273.0 and 273.5 K, indicating a melting point of ∼273.3 K for ice Ih. This predicted melting point is in excellent agreement with the experimental melting point of 273.15 K at ambient pressure. Also, the free energy curve shows a sharp change at the same temperature.

FIG. 9.

The (a) absolute entropy and (b) free energy of water as a function of temperature from the 2PT analysis of the MD trajectory. We find a sharp discontinuity between T = 273.0 and 273.5 K, corresponding to the melting point of ice Ih.

FIG. 9.

The (a) absolute entropy and (b) free energy of water as a function of temperature from the 2PT analysis of the MD trajectory. We find a sharp discontinuity between T = 273.0 and 273.5 K, corresponding to the melting point of ice Ih.

Close modal

We used two methods to evaluate the liquid density of water at 298 K and 1 atm. First, we performed MD simulations in the NPT ensemble for 1 nanosecond (ns). Since large fluctuations of the pressure during MD-NPT simulations can sometimes skew the NPT ensemble averages, we also computed the EOS of the liquid water at 298 K for densities in the range of 0.98 to 1.02 gr/cm3. For each density, we utilized MD-NVT simulations at 298 K for 1 ns and averaged the pressure values over the last 0.5 ns. The pressure versus density plot is shown in Fig. 10, leading to a density of 0.9965 gr/cm3 at 1 atm in perfect agreement with the experimental density of 0.9965 gr/cm3. Our MD-NPT also resulted in a density of 0.9960 gr/cm3 also in excellent agreement with experimental density (see Fig. S8 of the supplementary material).

FIG. 10.

Equation of state of water predicted at 298 K. The data are fitted to a second order polynomial (dashed line), leading to a density of water at 1 atm of 0.9965 gr/cm3 at 298 K, in perfect agreement with the experimental density of 0.9965 gr/cm3.

FIG. 10.

Equation of state of water predicted at 298 K. The data are fitted to a second order polynomial (dashed line), leading to a density of water at 1 atm of 0.9965 gr/cm3 at 298 K, in perfect agreement with the experimental density of 0.9965 gr/cm3.

Close modal

We computed the enthalpy of vaporization (ΔHvap) as the difference between the gas (g) and liquid (l) enthalpies,

ΔHvapT=Hp,TgHp,Tl,
(30)

which can be written in terms of gas and liquid total energies and volumes (Vg and Vl),

ΔHvapT=UgTUlT+pVgVl,
(31)

where Ug(T) and Ul(T) are the total energies of the gas and liquid phase, respectively. Equation (31) is further simplified with the assumptions that (i) VlVg, (ii) gas phase is ideal (i.e., pVg=RT), and (iii) the kinetic energy of water molecules in the liquid and gas phases is equal,

ΔHvapT=EgTElT+RT,
(32)

where Eg(T) and El(T) are the potential energies of the gas and liquid phase, respectively. It should be noted that Eq. (32) requires no correction for polarization effects since the PQEq model explicitly includes the polarization energy. Using Eq. (32), we computed ΔHvap = 10.36 kcal/mol at 298 K compared to be in excellent agreement with experimental ΔHvap = 10.52 kcal/mol.

We computed the S0 of liquid at 298 K using 2PT methodology. For this calculation, we extracted the last 100 ps trajectory of the 1 ns MD-NVT simulation described above. The 2PT analysis was performed for five different 20 ps intervals and averaged over the 5 cycles. This leads to S0 = 68.43 (J/mol)/K, just 2% below that the experimental value of S0 = 69.9 (J/mol)/K.68 This accurate value from RexPoN for such absolute thermodynamics properties as entropy suggests that this model will reproduce a wide variety of properties of water. The computed S0 is usually underestimated by other water models. Thus AIMD-PBE (−27%), ReaxFF (−13%), TIP4P-2005 (−18%), TIP4P-ice (−25%), SPC/Fw (−13%), F3C (−12%), effective FF from AIMD (−6%), and TIP4/F (−16%) are compared to experiment.64 

We computed the static dielectric constant (ε) from dipole moment fluctuations

ε=ε+4π3kBTVM2M2,
(33)

where kB is the Boltzmann constant, T is the temperature, V is the volume, M is the total dipole moment of the simulation cell, and ε is the high frequency or optical dielectric constant. The PQEq formalism allows the direct calculation of ε from the classical fluctuation of the electrons (shells) as a measure of the total polarizability of the simulation cell.17 It can be calculated from the dipole moment fluctuations due to the shell movement around the cores using

ε=1+4π3kBTVMs2Ms2,
(34)

where Ms indicates the dipole moment resulted from the shell displacement from its core (Ms=iNΔri). To compute the dielectric constant, we performed = 1 ns MD-NVT simulations at 298 K, leading to ε0 = 76.1 and ε = 1.31, in good agreement with the experimental dielectric constant69 value of 78.4 and the optical dielectric constant70 of 1.79. The change with time of the ε, M2, and M2 are shown in Fig. 11. As shown in this figure, after 300 ps, the values of ε and M2 reach to a plateau and M2 becomes negligible.

FIG. 11.

The change of static dielectric constant (ε) and dipole moment fluctuations (M2 and M2) with time during 1 ns MD-NVT simulations at T = 298 K and ambient pressure. After 300 ps, the values of ε and M2 reach to a plateau and M2 becomes negligible.

FIG. 11.

The change of static dielectric constant (ε) and dipole moment fluctuations (M2 and M2) with time during 1 ns MD-NVT simulations at T = 298 K and ambient pressure. After 300 ps, the values of ε and M2 reach to a plateau and M2 becomes negligible.

Close modal

The molecular distances and orientations in water are mainly determined by the HB network within the first shell of each molecule, making the structural analysis of the first nearest neighbors crucial. We computed the oxygen-oxygen radial distribution function (gOO) and coordination number (NOO) of water at 298 K and ambient pressure over 1 ns MD-NVT simulations. The results together with the experimental data are shown in Fig. 12. The experimental data are based on neutron scattering (Exp1) and combination of X-ray and Neutron experiments (Exp2).71,72 The first peak of RexPoN (2.84 Å) is in excellent agreement with Exp1 (2.85 Å) and Exp2 (2.76 Å). Most empirical water models were adjusted to have a first peak at ∼2.76 Å, but they lead to a sharp first peak.72–74 Figure 12 includes the RDF for the TIP4P/2005 FF74 (a typical empirical FF) and several ab initio based models including SCAN,75 CC-pol, and MB-pol potentials. All these models lead to a poor fit to the experimental gOO. Although the ab initio models lead to a much better gOO than the empirical FF. By contrast, RexPoN gives a slope and height of the first peak in excellent agreement with experimental data. Our computed second peak (4.2 Å) is also within the experimental range of 4.2 to 4.7 Å. As shown in Fig. 12, RexPoN leads to a first minimum of the goo at 3.43 Å, leading to Noo = 4.7, in excellent agreement with experiment. This value shows more than four water molecules in the first shell of water, suggesting a broken tetrahedral structure for water.

FIG. 12.

Comparison of the oxygen-oxygen radial distribution function at T = 298 K and ambient pressure between experiment, RexPoN, TIP4P-2005f,74 SCAN,78 CC-pol,20 and MB-pol21 water models. The experimental data are based on neutron scattering (Exp1)71 and combination of X-ray and neutron experiments (Exp2).72 

FIG. 12.

Comparison of the oxygen-oxygen radial distribution function at T = 298 K and ambient pressure between experiment, RexPoN, TIP4P-2005f,74 SCAN,78 CC-pol,20 and MB-pol21 water models. The experimental data are based on neutron scattering (Exp1)71 and combination of X-ray and neutron experiments (Exp2).72 

Close modal

The self-diffusion coefficient (Ds) was derived from the ionic mean-squared displacement (MSD) curve using the 3D diffusion equation of Einstein

MSDt=limtritri(t0)2=6Dst,
(35)

where ri(t) is the position of particle i at time t and t0 is the time origin and could be any time during the simulation. To provide better sampling, we computed an averaged MSD at time t with respect to all of the frames smaller than t. The values of MSD were computed every 1 ps during 1 ns NVT simulations at 298 K and ambient pressure. The change of MSD with time along with its log-log curve is shown in Fig. 13. The log-log curve shows that by 300 ps, the system has reached the Fickian regime with a slope of 1.0, satisfying logMSD=logt+log6Ds. Our computed Ds = 4.21 × 10−5 cm2 s−1 is higher than the experimental value76 of 2.27 × 10−5 cm2 s−1. From 2PT analysis, we also computed the diffusivity using the density of states at zero frequency, DoS(0) = 12 mN DS/(kT), which leads to Ds = 4.23 × 10−5 cm2 s−1. This differs the MSD value by only 0.5%, suggesting that the 2PT analysis can be used to compute the diffusivity of large systems from short simulation times. The larger computed value of diffusivity by RexPoN is puzzling to us. We have done this both for long time scales and from 2PT and obtained the same results. This discrepancy is puzzling because RexPoN leads to good accuracy for the entropy of water [S0 = 68.43 (J/mol)/K compared to the experimental value of 69.90 (J/mol)/K]. It is known that diffusivity and entropy of a system are directly related to each other. Therefore, one would expect RexPoN to give an accurate value for the diffusivity as it describes the entropy of the system very well.

FIG. 13.

The change of mean-squared displacement with time during 1 ns MD-NVT simulations at 298 K and 1 atm pressure. The slope of the fitted red line is used to determine the self-diffusivity coefficient (Ds) of the liquid according to the Einstein equation. The slope of the log-log curve (inset) becomes equal to 1.0 after about 300 ps.

FIG. 13.

The change of mean-squared displacement with time during 1 ns MD-NVT simulations at 298 K and 1 atm pressure. The slope of the fitted red line is used to determine the self-diffusivity coefficient (Ds) of the liquid according to the Einstein equation. The slope of the log-log curve (inset) becomes equal to 1.0 after about 300 ps.

Close modal

The RexPoN force field based entirely on high quality QM calculations leads to excellent agreement with the experimental properties of the solid and liquid phases of water. This attests to the quality of the QM calculations and to the validity of the RexPoN formalism. We are hopeful that the extension of the RexPoN methodology to other systems, fitting the parameters from QM, may lead to a new generation of FF for very accurate predictions on electrolytes and other materials. We believe that the RexPoN FF for water could provide a significant reduction in the costs for full solvent electrochemical QM calculations by allowing the QM water to be restricted to just a couple of layers.

The supplementary material contains the comparison between PQEq and QM charge distributions in the water dimer, the crystal structures of solid hydrogen and solid oxygen, two-body vdW energy curves, energy changes with the volume of the hydrogen crystal, bond energy to the bond order relationship, structures of 10 low-lying stationary points of the water dimer, global minimized structures of 19 water clusters, comparison between RexPoN and QM for the 4 scan cases, comparison between RexPoN and QM energies for water dimers and clusters, the liquid density evaluation using MD-NPT simulations, and the computational cost of RexPoN.

This work was supported by the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC00014607. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation Grant No. ACI-1548562. W.A.G. was supported by NSF (CBET 1512759).

1.
J. F.
Ouyang
and
R.
Bettens
,
CHIMIA Int. J. Chem.
69
(
3
),
104
111
(
2015
).
2.
B.
Guillot
,
J. Mol. Liq.
101
(
1-3
),
219
260
(
2002
).
3.
F.
Prielmeier
,
E.
Lang
,
R.
Speedy
, and
H.-D.
Lüdemann
,
Phys. Rev. Lett.
59
(
10
),
1128
(
1987
).
4.
P. H.
Poole
,
F.
Sciortino
,
U.
Essmann
, and
H. E.
Stanley
,
Nature
360
(
6402
),
324
(
1992
).
5.
C. G.
Salzmann
,
P. G.
Radaelli
,
A.
Hallbrucker
,
E.
Mayer
, and
J. L.
Finney
,
Science
311
(
5768
),
1758
1761
(
2006
).
6.
C. G.
Salzmann
,
P. G.
Radaelli
,
E.
Mayer
, and
J. L.
Finney
,
Phys. Rev. Lett.
103
(
10
),
105701
(
2009
).
7.
W. L.
Jorgensen
,
J. Am. Chem. Soc.
103
(
2
),
335
340
(
1981
).
8.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
,
J. Chem. Phys.
79
(
2
),
926
935
(
1983
).
9.
M. W.
Mahoney
and
W. L.
Jorgensen
,
J. Chem. Phys.
112
(
20
),
8910
8922
(
2000
).
10.
H. J.
Berendsen
,
J. P.
Postma
,
W. F.
van Gunsteren
, and
J.
Hermans
,
Intermolecular Forces
(
Springer
,
1981
), pp.
331
342
.
11.
M. W.
Mahoney
and
W. L.
Jorgensen
,
J. Chem. Phys.
115
(
23
),
10758
10768
(
2001
).
12.
Y.
Wu
,
H. L.
Tepper
, and
G. A.
Voth
,
J. Chem. Phys.
124
(
2
),
024503
(
2006
).
13.
S. W.
Rick
,
S. J.
Stuart
, and
B. J.
Berne
,
J. Chem. Phys.
101
(
7
),
6141
6156
(
1994
).
14.
C.
Millot
and
A. J.
Stone
,
Mol. Phys.
77
(
3
),
439
462
(
1992
).
15.
E. M.
Mas
,
K.
Szalewicz
,
R.
Bukowski
, and
B.
Jeziorski
,
J. Chem. Phys.
107
(
11
),
4207
4218
(
1997
).
16.
C. J.
Burnham
,
J.
Li
,
S. S.
Xantheas
, and
M.
Leslie
,
J. Chem. Phys.
110
(
9
),
4566
4581
(
1999
).
17.
G.
Lamoureux
,
A. D.
MacKerell
, Jr.
, and
B.
Roux
,
J. Chem. Phys.
119
(
10
),
5185
5197
(
2003
).
18.
B. T.
Thole
,
Chem. Phys.
59
(
3
),
341
350
(
1981
).
19.
X.
Huang
,
B. J.
Braams
, and
J. M.
Bowman
,
J. Phys. Chem. A
110
(
2
),
445
451
(
2006
).
20.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
van der Avoird
,
J. Chem. Phys.
128
(
9
),
094313
(
2008
).
21.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
(
12
),
5395
5403
(
2013
).
22.
M.
Riera
,
N.
Mardirossian
,
P.
Bajaj
,
A. W.
Götz
, and
F.
Paesani
,
J. Chem. Phys.
147
(
16
),
161715
(
2017
).
23.
M.
Sprik
,
J.
Hutter
, and
M.
Parrinello
,
J. Chem. Phys.
105
(
3
),
1142
1152
(
1996
).
24.
C. M.
Handley
,
G. I.
Hawe
,
D. B.
Kell
, and
P. L.
Popelier
,
Phys. Chem. Chem. Phys.
11
(
30
),
6365
6376
(
2009
).
25.
A.
Shank
,
Y.
Wang
,
A.
Kaledin
,
B. J.
Braams
, and
J. M.
Bowman
,
J. Chem. Phys.
130
(
14
),
144314
(
2009
).
26.
M.
Ceriotti
,
W.
Fang
,
P. G.
Kusalik
,
R. H.
McKenzie
,
A.
Michaelides
,
M. A.
Morales
, and
T. E.
Markland
,
Chem. Rev.
116
(
13
),
7529
7550
(
2016
).
27.
S.
Habershon
,
T. E.
Markland
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
131
(
2
),
024501
(
2009
).
28.
S.
Naserifar
,
D. J.
Brooks
,
W. A.
Goddard
 III
, and
V.
Cvicek
,
J. Chem. Phys.
146
(
12
),
124117
(
2017
).
29.
J. J.
Oppenheim
,
S.
Naserifar
, and
W. A.
Goddard
 III
,
J. Phys. Chem. A
122
(
2
),
639
645
(
2018
).
30.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
(
15
),
154104
(
2010
).
31.
S.
Grimme
,
S.
Ehrlich
, and
L.
Goerigk
,
J. Comput. Chem.
32
(
7
),
1456
1465
(
2011
).
32.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
(
18
),
3865
(
1996
).
33.
A. D.
Becke
,
J. Chem. Phys.
98
(
7
),
5648
5652
(
1993
).
34.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
(
2
),
785
(
1988
).
35.
C.
Barrett
,
L.
Meyer
, and
J.
Wasserman
,
J. Chem. Phys.
47
(
2
),
592
597
(
1967
).
36.
Y. A.
Freiman
and
H.-J.
Jodl
,
Phys. Rep.
401
(
1-4
),
1
228
(
2004
).
37.
P.
Loubeyre
,
R.
LeToullec
,
D.
Hausermann
,
M.
Hanfland
,
R.
Hemley
,
H.
Mao
, and
L.
Finger
,
Nature
383
(
6602
),
702
(
1996
).
38.
J.
Kohanoff
,
S.
Scandolo
,
G. L.
Chiarotti
, and
E.
Tosatti
,
Phys. Rev. Lett.
78
(
14
),
2783
(
1997
).
39.
M.
Städele
and
R. M.
Martin
,
Phys. Rev. Lett.
84
(
26
),
6070
(
2000
).
40.
41.
L. F.
Lundegaard
,
G.
Weck
,
M. I.
McMahon
,
S.
Desgreniers
, and
P.
Loubeyre
,
Nature
443
(
7108
),
201
(
2006
).
42.
G.
Weck
,
S.
Desgreniers
,
P.
Loubeyre
, and
M.
Mezouar
,
Phys. Rev. Lett.
102
(
25
),
255503
(
2009
).
43.
Y.
Akahama
,
H.
Kawamura
, and
O.
Shimomura
,
Phys. Rev. B
64
(
5
),
054105
(
2001
).
44.
R.
Dovesi
,
R.
Orlando
,
A.
Erba
,
C. M.
Zicovich-Wilson
,
B.
Civalleri
,
S.
Casassa
,
L.
Maschio
,
M.
Ferrabone
,
M.
De La Pierre
, and
P.
D’Arco
,
Int. J. Quantum Chem.
114
(
19
),
1287
1317
(
2014
).
45.
J.
Heyd
,
J. E.
Peralta
,
G. E.
Scuseria
, and
R. L.
Martin
,
J. Chem. Phys.
123
(
17
),
174101
(
2005
).
46.
H.
Kim
,
J.-M.
Choi
, and
W. A.
Goddard
 III
,
J. Phys. Chem. Lett.
3
(
3
),
360
363
(
2012
).
47.
Y.
Liu
and
W. A. I.
Goddard
,
Mater. Trans.
50
(
7
),
1664
1670
(
2009
).
48.
L.
Liu
,
Y.
Liu
,
S. V.
Zybin
,
H.
Sun
, and
W. A.
Goddard
 III
,
J. Phys. Chem. A
115
(
40
),
11016
11022
(
2011
).
49.
A. K.
Rappé
,
C. J.
Casewit
,
K.
Colwell
,
W.
Goddard
 III
, and
W.
Skiff
,
J. Am. Chem. Soc.
114
(
25
),
10024
10035
(
1992
).
50.
51.
A. C.
Van Duin
,
S.
Dasgupta
,
F.
Lorant
, and
W. A.
Goddard
,
J. Phys. Chem. A
105
(
41
),
9396
9409
(
2001
).
52.
J. S.
Sims
and
S. A.
Hagstrom
,
J. Chem. Phys.
124
(
9
),
094101
(
2006
).
53.
L.
Bytautas
and
K.
Ruedenberg
,
J. Chem. Phys.
132
(
7
),
074109
(
2010
).
54.
F. R.
Gilmore
,
J. Quant. Spectrosc. Radiat. Transfer
5
(
2
),
369
390
(
1965
).
55.
A.
Kumar
,
S. R.
Gadre
,
N.
Mohan
, and
C. H.
Suresh
,
J. Phys. Chem. A
118
(
2
),
526
532
(
2014
).
56.
J. T.
Su
,
X.
Xu
, and
W. A.
Goddard
,
J. Phys. Chem. A
108
(
47
),
10518
10526
(
2004
).
57.
Schrödinger Release 2015-2: Jaguar, version 8.8, Schrödinger, LLC, New York, NY,
2015
.
58.
F. H.
Stillinger
,
Science
209
(
4455
),
451
457
(
1980
).
59.
S.
Plimpton
,
J. Comput. Phys.
117
(
1
),
1
19
(
1995
).
60.
W. G.
Hoover
,
Phys. Rev. A
31
(
3
),
1695
(
1985
).
61.
S.
Nosé
,
J. Chem. Phys.
81
(
1
),
511
519
(
1984
).
62.
J.-P.
Ryckaert
,
G.
Ciccotti
, and
H. J.
Berendsen
,
J. Comput. Phys.
23
(
3
),
327
341
(
1977
).
63.
S. K.
Reddy
,
S. C.
Straight
,
P.
Bajaj
,
C.
Huy Pham
,
M.
Riera
,
D. R.
Moberg
,
M. A.
Morales
,
C.
Knight
,
A. W.
Götz
, and
F.
Paesani
,
J. Chem. Phys.
145
(
19
),
194504
(
2016
).
64.
T. A.
Pascal
,
D.
Schärf
,
Y.
Jung
, and
T. D.
Kühne
,
J. Chem. Phys.
137
(
24
),
244507
(
2012
).
65.
C.
Vega
and
J. L.
Abascal
,
Phys. Chem. Chem. Phys.
13
(
44
),
19663
19688
(
2011
).
66.
S.-T.
Lin
,
P. K.
Maiti
, and
W. A.
Goddard
 III
,
J. Phys. Chem. B
114
(
24
),
8191
8198
(
2010
).
67.
T. A.
Pascal
,
S.-T.
Lin
, and
W. A.
Goddard
 III
,
Phys. Chem. Chem. Phys.
13
(
1
),
169
181
(
2011
).
68.
J.
Cox
,
D.
Wagman
, and
V.
Medvedev
, USA Google Scholar, New York, 1989.
69.
D. P.
Fernández
,
Y.
Mulev
,
A.
Goodwin
, and
J. L.
Sengers
,
J. Phys. Chem. Ref. Data
24
(
1
),
33
70
(
1995
).
70.
A.
Buckingham
,
Proc. R. Soc. A
238
(
1213
),
235
244
(
1956
).
71.
A.
Narten
,
W.
Thiessen
, and
L.
Blum
,
Science
217
(
4564
),
1033
1034
(
1982
).
72.
A.
Soper
,
J. Phys.: Condens. Matter
19
(
33
),
335206
(
2007
).
73.
H.
Berendsen
,
J.
Grigera
, and
T.
Straatsma
,
J. Phys. Chem.
91
(
24
),
6269
6271
(
1987
).
74.
M. A.
González
and
J. L.
Abascal
,
J. Chem. Phys.
135
(
22
),
224516
(
2011
).
75.
M.
Chen
,
H.-Y.
Ko
,
R. C.
Remsing
,
M. F. C.
Andrade
,
B.
Santra
,
Z.
Sun
,
A.
Selloni
,
R.
Car
,
M. L.
Klein
, and
J. P.
Perdew
,
Proc. Natl. Acad. Sci. U. S. A.
114
(
41
),
10846
10851
(
2017
).
76.
R.
Mills
,
J. Phys. Chem.
77
(
5
),
685
688
(
1973
).
77.
C.
McBride
,
C.
Vega
,
E. G.
Noya
,
R.
Ramírez
, and
L. M.
Sesé
,
J. Chem. Phys.
131
(
2
),
024506
(
2009
).
78.
J.
Wiktor
,
F.
Ambrosio
, and
A.
Pasquarello
,
J. Chem. Phys.
147
,
216101
(
2017
).

Supplementary Material