Classical water potential models take up the majority and have played a dominant role in molecular dynamics simulations of water. Rigid models, e.g., TIP4P, TIP4P-Ew [5], TIP4P/ice [15], TIP4P/2005 [6], TIP5P [23], TIP5P-Ew [8], and NvdE [7], have been studied in detail from various perspectives. Melting temperatures of these models obtained by different researchers are not exactly same but within the uncertainty. Hence, melting points of the mentioned models are adapted qualitatively from the published literature and the form “model-melting point (uncertainty)” is used for a simple expression. Melting temperatures are listed as TIP4P-230.5 (3) K, TIP4P-Ew-244 (3) K, TIP4P/ice-270 (3) K, TIP4P/2005-250.5 (3) K, TIP5P-272 (3) K, TIP5P-Ew-271 (3) K [24], and NvdE-289 (3) K [21] (see Figure 1). For complete comparisons with bulk water, the melting points of 20-molecule water cluster using TIP4P, TIP4P/ice and TIP4P/2005 were estimated within 150 ± 20 K. TIP4P/ice had the highest value followed by TIP4P/2005 and TIP4P [18], which is consistent with bulk water. Deviations from experiments are justified by distinguishing the position of negative charges and the resulting quadrupole moment of these models. One negative charge is located at the bisector of the two OH bonds in TIP4P-like models, while two negative charges in TIP5P and TIP5P-Ew are located at the lone-pair electron sites of the oxygen atom. As for NvdE model, three negative charges are located at the lone-pair sites and the bisector simultaneously. In three-charge models, once the quadrupole moment is closer to the real value, the melting point is more consistent with the experimental value [25]. However, the melting temperature of four-charge models is correlated with a stable tetrahedral distribution of charges. This type of charge distribution enhances the hydrogen bond network providing a higher melting point for the NvdE model compared to TIP4P model. Accounting for the geometry of NvdE which combines two methodologies, the negative charge distribution is inferred to dramatically influence the melting temperature. Conversely, parameters of rigid models are set to match experimental data of target properties at different conditions. For instance, TIP4P/ice was specifically designed to calculate solid phase properties [15], and TIP4P/2005 was designed to reproduce the temperature of maximum density and other common properties [6], although the two models have the similar geometry and charge distribution. Consequently, estimations of anticipated properties like the melting temperature can be guaranteed for specific models, but predictions of other properties mostly fail for exceeding the fit range. Finally, measurements of melting points are associated with phase transitions that are sensitive to pressure and temperature. As the pressure increases, low-density water transforms continuously to high-density water due to interstitial molecules, which was simulated by TIP4P [26]. The relation that a higher pressure corresponds to a lower melting point was proved and the validity of TIP4P was checked concomitantly.

The validity of rigid models is unsurprisingly questionable since they ignore flexible and polarizable features. The flexibility affects intramolecular interactions in terms of variable OH bond lengths and angles. The flexible q-TIP4P/F [27] model with OH stretches described by Morse-type functions determined the melting temperature to be 259 (1) K. The inclusion of such intramolecular interaction indicated no obvious improvements to the melting temperature. In addition to flexibility, the polarizability can reveal the potential to respond instantaneously to heterogenous environments and charge transformations. The TIP4P-FQ model [28] which modifies TIP4P through fluctuating charges contains a self-polarization term, causing stronger interactions between molecules. Then, the corresponding melting temperature was found to be 303 (8) K [29]. Another polarizable model named SWM4-NDP [30] applies the Drude oscillator technique to represent polarizability and has the quadrupole moment assigned at approximately the experimental value. Surprisingly, the melting temperature of SWM4-NDP was 185 (10) K, which was much lower than 273.15 K [31]. The polarizable POL3 [32] model yielded a melting point of 180 (10) K [17] and thus its suitability for ice simulation is limited. In this model, induced dipoles are determined by point dipoles located at the three atom sites and the surrounding electric field. Although induced point dipoles increase the total dipole moment, the total quadrupole moment is rather small. The POL4D model [9] with Rowlinson five-site geometry is polarized in the simplest isotropic linear way. The melting temperature of this model was estimated as 260 (2) K. Importantly, methods treating polarizability and structures of polarizable models can affect the simulated melting temperature, which is clearly demonstrated by deviations from the experimental value [31] (see Figure 1). In most cases, the addition of polarizability tends to lower the melting temperature with an exception of TIP4P-FQ. When the polarizability is introduced based on common rigid models, this modification may destabilize the solid structure and increase the disorder of molecules, which accounts for the lower melting temperature.

In contrast to roughness of empirical potential models, ab initio potential models express the sophisticated representation of Born-Oppenheimer potential energy surface. TTM2.1-F [11] which employed gas phase monomer potential energy and dipole moment surfaces previously yielded a melting value of 242.5 (1.5) K [33]. Later, the flexible and polarizable TTM3-F model [16] having similar structural and dynamical characteristics with the previous version converged at 248 (2.5) K [34]. The highlight of the latter model is the ability to reproduce vibrational spectra both in water clusters and liquid water, which are congruent with experiments [16]. To some extent, this agreement ensures the validity of the designed model. The EFP model [35] is derived from ab initio quantum chemistry with no empirical parameters and water molecules are represented as combinations of five interaction energy terms, including the Coulombic, polarization, exchange repulsion, dispersion and charge transfer interactions. The estimated melting point was 381 (15) K [20] which was substantially higher than the experimental value. The iAMOEBA model [36] is parameterized using the ForceBalance method which utilizes experimental data and high-level ab initio calculation results. Electronic polarizability is treated in a direct approximation, which reduces the computational cost and maintains an accurate description. The melting point was calculated to be 261 (2) K, not far from the experimental value. The WAIL model was created as a problem-specific force field for water using the adaptive force matching method with only electronic structure information as input. Thus, this model determined a melting point of 270 (2.5) K [37]. Although ab initio potential models are desired to predict macroscopic properties and interactions well, the availability of such models that can predict melting temperature is still problematic (see Figure 1). However, information derived from the electronic structure and potential energy surface is encouraging for a deeper understanding of water molecules. Despite this, how to use such information of intermolecular interactions to account for the melting properties remains challenging, which is clearly shown by the rather high melting temperature of the EFP model. TTM2.1-F and TTM3-F that rely on water clusters interactions may underestimate bulk interactions and lead to lower melting temperatures, while iAMOEBA which employs hybrid data and WAIL which uses a quantum mechanics/molecular mechanics procedure, perform well about the melting prediction.

Great caution should be exercised when DFT-based molecular dynamics simulations are performed. Most exchange-correlation functionals fail to predict physical and thermodynamic properties of water, resulting from the lack of delicate interactions. In view of the complexity and accuracy of functionals, Perdew-Burke-Ernzerhof (PBE) and Becke-Lee-Yang-Parr (BLYP) were chosen to determine the melting temperature via Born-Oppenheimer molecular dynamics simulations [14]. When the prerequisite density (ρ = 1 g/cm³) was satisfied, the average pressure of P = 2,500 bar (PBE) and P = 10,000 bar (BLYP) were much higher than the ambient pressure. Melting temperatures corresponding to these unusual pressures were 417(3) K (PBE) and 411 (4) K (BLYP). Supercooled water and overstructured liquid were predicted at ambient conditions, which strongly suggests that pure functionals cannot effectively describe accurate interactions. To this end, an empirical van der Waals (vdW) potential correction for long-range dispersion effects has been verified to dramatically improve descriptions of molecular complexes [38]. The dispersion correction plays a significant role in improving simulation results, yielding superior water densities, radial distribution functions and self-diffusion coefficients. BLYP-D functional [34] estimated the melting temperature as 360 (10) K, showing noticeable progress compared to the original functional. The development of approximate DFT approaches is an active area of research for the chemical accuracy in large systems, such as dispersion-corrected atom-centered potential, nonlocal van der Waals functional and pure density functional [39]. However, BLYP functional with D3 corrections provided a melting value of 325 K [22] which is still higher than the experimental value. Besides, effects of vdW interactions on properties of water are explored by using neural network potentials which produce the accurate potential energy surface. Calculated melting temperatures of revised Perdew–Burke–Ernzerhof (RPBE) and BLYP with and without D3 corrections are displayed as BLYP-323 (3) K, BLYP-vdW-283 (2) K, RPBE-267 (2) K and RPBE-vdW-274 (3) K [40]. The inclusion of vdW forces not only resulted in more realistic structures but also produced a shift of melting temperature toward 273.15 K. Notably, a discrepancy occurs when values of BLYP are compared with earlier reports in term of different settings, but the diminution of the melting point due to vdW corrections is consistent, which further proves the role of vdW corrections. Clearly, DFT-based simulations roughly estimate the melting point because of the inherent complexity and shortcomings of functionals. Higher functionals like hybrid functionals or functionals with advanced dispersion corrections are suggested for excellent performance.

Monte Carlo simulations with second-order Møller-Plesset (MP2) perturbation theory have been performed for understanding structures and properties of liquid water [41]. Similar to DFT-based simulations, MP2-based simulations are also a quantum mechanical method. Moreover, MP2 theory accounts for dispersion interactions and recovers many dynamic correlations. Simulations for the melting temperature of water with MP2 are still missing, but the density of liquid water at ambient conditions (T = 295 K and P = 1 bar) has been obtained, i.e., 1.02 g/mL [41]. To obtain the correct density, DFT with a generalized gradient approximation needs a higher temperature and pressure than MP2. Polarizability and dispersion interactions were found to explain these differences. The TTM3-F model with dispersion and polarizability scale factors was chosen to perform simulations individually to test effects of the dispersion and polarizability [42]. A smaller dispersion interaction leads to less dense liquid water, which indicates that molecules would move faster. Then, a larger pressure is needed on the system to maintain the density at 1 g/cm³. In fact, the correction of vdW interactions in BLYP potential eliminates the demand for higher pressures [34]. Furthermore, a larger polarizability was found to result in overstructured water, which would melt at a higher temperature. Enhanced tetrahedral hydrogen bonds prevent the motion and rotation of molecules at ambient conditions and more energy is needed to break these restrictions [42]. Through the comparison of dispersion and polarizability, the melting point is further explained by the effect of the intermolecular forces which are determined by the essential description of the molecular model.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.