Water is commonly referred to as the universal solvent because it dissolves more substances than any other known liquid. It is also essential for life and understanding and simulating its role in all biological processes is extremely important. The exceptional solvating capacity of water is thought to result from its unique molecular properties, including high dielectric constant, hydrogen bonding structures, mobility, self-ionization, and strong permanent dipole moment. The combination of these unique properties allows water to facilitate many solvation and self-assembly phenomena. As such, it is a critical component of molecular dynamics (MD) simulations of proteins, 2-5 DNA, membranes, polymers, in solvent mixtures involving water, and water-solvated materials. MD simulations are a powerful tool for providing atomistic insight into experimental systems of interest. Moreover, these techniques are accessible to researchers with little prior experience in modeling or computational approaches due to the development of freely available force field parameterization tools such as CGenFF, LigParGen, and others, as well as tools for constructing systems and input files such as VMD and CHARMM-GUI. While this is undeniably of great benefit to the scientific community at large, the number and diversity of water models currently available (Figure ) means it is critical that new users are aware of the potential pitfalls and concerns regarding the choice of water model for their MD simulations. For example, reproducing the correct dielectric constant of the individual components is essential for correctly describing the properties of a solvent mixture, and reproducing the surface tension of water is necessary for describing surfactants, monolayers, and interfaces. For transport properties such as hydrogen bond lifetimes and viscosity, a reasonable value of the self-diffusion coefficient of water is desirable. To ensure that simulated systems involving water are reflective of real-world phenomena, it is crucial to validate that the chosen model can reproduce the intrinsic properties of water. Here, we use the term 'water model' to describe the complete set of force field parameters (i.e. bond lengths, angles, force constants for flexible models, atomic partial charges, and Lennard-Jones sigma and epsilon values that describe van der Waals interactions) necessary to perform MD simulations of water. While simulations of water have been reported since the 1960s, the first generation of water models that are still commonly in use today was developed in the 1980s. The Simple Point Charge (SPC) water model, with a geometry based on the tetrahedral structure of water in ice and published by Berendsen et al. in 1981, is typically used with the GROningen MOlecular Simulation (GROMOS) force field. The Transferable Interaction Potential 3-point model (TIP3P), with a geometry based on gas-phase water and published by is still used with both the Assisted Model Building with Energy Refinement (AMBER) and Optimized Potentials for Liquid Simulations (OPLS) force fields. The slight differences in bond lengths and angles allow for SPC-based models to reproduce the experimentally observed second peak in the O-O radial distribution function, whereas the second peak is missing entirely from the initial TIP3P model. Radial distribution functions, which reflect molecular structure in terms of density, are commonly used for comparison as they can be directly compared to experimental data extracted from X-ray scattering of molecular fluids. Further modifications to the TIP3P water model, including the addition of Lennard-Jones potentials to the hydrogens, resulted in the CHARMM TIP3P water model (also referred to as mTIP3P or TIPS3P, and referred to as TIPS3P hereafter in this work) which is used with the Chemistry at HARvard Macromolecular Mechanics (CHARMM) force field. In the same 1983 paper as the original TIP3P model, Jorgensen et al. also described a 4-point model, TIP4P, where the negative charge of the oxygen was offset from the Lennard-Jones potential, resulting in improved structural properties. The TIP4P water model is also commonly used with the OPLS force field. For the Consistent-Valence Force Field (CVFF), the CVFF water model has a geometry similar to TIP3P and non-bond parameters similar to SPC. It was developed concurrently with the rest of the force field, and is notable for having a flexible geometry, the use of a Morse potential for the O-H bond, and cross-terms to reproduce vibrational frequencies. These models remain a popular choice due to their simple geometry, low computational cost, and relative accuracy. With the exception of CVFF, these models have a rigid geometry, which enables the use of a larger (e.g. 2 fs) timestep in MD simulations. While these models reproduced some experimental properties of water reasonably well, limited computational resources available at the time necessitated small system sizes of 125 to 216 water molecules, short simulation times of less than 20 ps, and short cut-offs of non-bond interactions typically between 7.5 and 8.5 Å. Moreover, many structural properties, such as densities and dielectric constants, were calculated using Monte Carlo, rather than MD, simulations. Monte Carlo simulations rely on statistical mechanics rather than dynamics, and may not reproduce the collective motion of molecules the same way as MD simulations. In the 1990s, Class II force fields based on ab initio potential energy surfaces were developed. Like CVFF, they include cross terms to better reproduce experimental vibrational frequencies, and in contrast to Class I force fields, use a 9-6 rather than 12-6 functional form for the Lennard-Jones interactions (Equations S1 and S2, Supporting Information). The Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) and Polymer Consistent Force Field (PCFF) force fields are the two Class II force fields most commonly still in use, primarily for polymers and materials applications. Like CVFF, these force fields were developed to be used with commercial software, but the force field parameters for the COMPASS and PCFF water models were published in 2012 as part of the INTERFACE force field. While many of these water models performed well at 25 °C, few could accurately reproduce the temperature of maximum density of water of 4 °C. In 2000, Mahoney et al. published the 5-point TIP5P water model, which reproduced the density of water in the temperature range of -37 to 100 °C and a pressure range of 1 to 10,000 atm. Following the initial publication of these 3-point, 4-point, and 5-point models, many modifications have been made to improve the properties based on new computational techniques and/or improvements in computing power. One of the first modifications to the initial model resulted in the SPC/E water model, which improved the reproduction of intermolecular potential energy. In order to increase compatibility with Ewald summation techniques for electrostatic interactions, the TIP3P-Ew, 61 TIP4P-Ew, and TIP5P-Ew 63 models were developed in 2004. In 2005, TIP4P was reparametrized to reproduce experimental properties with MD settings more in line with those currently used today, resulting in the TIP4P/2005 model. Flexible versions of the rigid models were also created, including SPC/Fw, 65 TIP3P/Fw 66 and TIP4P/2005f. In 2006, the polarizable Simple Water Model with 4 sites and Negative Drude Polarizability (SWM4-NDP) was published. Other models were produced with the goal of maximizing a single property. The TIP4P/Ice model was developed to study the properties of ice, while the TIP4P/ε, 70 SPC/ε, 71 H2O-DC, and FBA/ε 73 models were all developed to reproduce the experimental dielectric constant. The SPC/Eb model was developed to improve translational and rotational diffusion. The TIP4P-D model was parameterized to better reproduce London dispersion interactions, and subsequently modified for the a99SBdisp force field for protein simulations. In the last decade, further additional generalized water models have been produced. The TIP3P-FB and TIP4P-FB models, based on the ForceBalance approach, were published in 2014, and updated as the TIP3P-ST and TIP4P-ST models to include surface tension in the optimization procedure in 2019. In 2018 an update of TIP5P was published as TIP5P-2018. In 2014, the "optimal" point charge 4-point model (OPC) was published, followed by a 3-point variant (OPC3) in 2016, where it was suggested that the accuracy limit of a 3-point rigid non-polarizable model had been reached.