Densities and radial distribution functions (RDFs) were determined from the 20 ns NPT simulations. Following this, the simulation box length was set to the average value from the NPT simulations and additional NVT simulations were run in triplicate for 20 ns. These simulations were used for the calculation of coordination number, self-diffusion coefficient, and dielectric constant using PyLAT.py. For 5-point and polarizable models, PyLAT.py did not calculate the correct self-diffusion coefficient, and for these models the calculation was performed using the mean-squared displacement of oxygen atoms from LAMMPS and the Einstein relation. For the calculation of coordination number, it was necessary to output trajectories containing only atoms (i.e. no lone pairs) for the 5-point models. For the calculation of dielectric constant, it was necessary to run additional NVT simulations with explicit 4-point models, as the implicit 4-point model trajectories would not produce the correct dipole moments for 4-point models. PyLAT.py also failed to calculate the correct dielectric constant for the polarizable model, which we calculated with our in-house code. For the determination of surface tension, the box was extended to 120 Å in the z dimension and NVT simulations were run for each model, with the surface tension determined from the pressure tensors calculated in LAMMPS by the virial method. The solvation free energy of model solutes into each water model can be probed using Widom's test particle insertion method. The solute is inserted into an equilibrated configuration of the pure solvent at position xs. The change in the potential energy ΔU upon insertion, without allowing the system to relax, is then calculated. For MD simulations in the NpT ensemble, the solvation free energy, ΔGsol, equivalent to the excess chemical potential, is given by