Given the sheer number of available water models and reported properties under a variety of disparate conditions (e.g., number of molecules, simulation time, non-bond cutoffs, simulation method, etc), choosing an appropriate water model for a particular study can be a daunting task, particularly for those new to MD simulation techniques. It is therefore desirable to have a comparison of water models where the properties were determined under the same conditions for every model. Here, we report the density, coordination number, surface tension, static dielectric constant, self-diffusion coefficient, and free energy of hydrophobic solvation of 30 representative 3-point, 4-point, and 5-point water models, and one polarizable model (Table ) using the open-source MD software LAMMPS and settings commonly used for simulations of biomolecules and materials today. We chose the SWM4-NDP model as a representative polarizable model because tools such as CHARMM-GUI and polarizer.py allow for the creation of LAMMPS input files for this model quickly and easily, even for non-expert users. This work serves as an update to the 2002 water model comparison by Guillot, and encompasses models published both before and after. While we don't examine continuum solvent models (i.e. implicit solvent) or coarse-grained models, we point the reader to a comprehensive review of water models published in 2017 by Onufriev and Izadi for further reading. Finally, we used machine learning methods to explore the contributions of the parameters (model type, partial charge, Lennard-Jones parameters, bond length and angle) to the surface tension, dielectric constant, and self-diffusion coefficient. Machine learning can be used to predict a wide variety of material, chemical, and biological properties, with several commercial and non-commercial open-source platforms that can be used to develop machine learning algorithms such as Schrödinger, SYBYL, 99 TensorFlow (Google), and BioPPSy. 101 Methods Molecular dynamics simulations. All simulations were performed with the MD code Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) release 5Jun19. LAMMPS was chosen as it can run Class I, Class II, and polarizable models easily, with input files generated by CHARMM-GUI, fftool, or Visual Molecular Dynamics (VMD). Moreover, LAMMPS provides the option of running 4-point water models as explicit 3-point models with the offset partial charge of the oxygen calculated internally. This allows 4-point models to be used from structure files that contain only 3point water (referred to hereafter as implicit 4-point simulations). LAMMPS is compatible with the PyLAT.py analysis tool, which facilitates the direct calculation of self-diffusion coefficients, dielectric constants, and radial distribution functions, among other properties. Initial configurations of 2,000 water molecules in a cubic box with approximate dimensions of 40 × 40 × 40 Å 3 were constructed using PACKMOL version 18.169, 105 and converted to LAMMPS data files using TopoTools 106 release 1.7 in VMD version 1.9.3. All 3-point and implicit 4-point models used the same initial configuration, with a conjugate gradient energy minimization step in LAMMPS to optimize the initial geometry for each model, using force constants of 500 kcal/Å² for the O-H bond and 50 kcal/radian² for the H-O-H angle with the rigid models, while flexible models used their default parameters. Explicit 4-point models, which were necessary for the calculation of dielectric constants, were constructed individually using PACKMOL, as were the 5-point models. The polarizable SWM4-NDP model was constructed as an explicit 4-point model using PACKMOL, with the Drude oscillator added by the Python tool polarizer.py. The full set of force field parameters for the water models employed in this study are shown in Table and Table for 3-point and 4-point models, respectively, Table for 5-point and the polarizable model, and Table for the force field parameters for flexible models. a In 4-point models the oxygen center of mass carries a charge of 0, and the charge is offset at point M (Figure ). For the 3-point and implicit 4-point simulations, initial random velocities with a Gaussian distribution at 1 K were assigned to each atom, and 500 ps of MD simulation was performed in the NVT (constant volume, constant temperature) ensemble with the system heated linearly from 1 K to 298.15 K. This was followed by 500 ps of NVT simulation at 298.15 K, and 4 ns of NPT (constant pressure, constant temperature) simulation at 298.15 K and 1.0 atm. Production runs were performed for 20 ns, with simulation snapshots saved every 10 ps.