To examine the model performance on bulk water properties, we performed both classical MD and PIMD simulations at 295K in NVT ensemble and computed the radial distribution functions (RDFs). Both classical and path-integral molecular dynamics were conducted using 1018 water molecules in a cubic box, with the density set to be 0.03338 molecule per 3 Å . 129 The intermolecular potential was computed using our models and the intramolecular part was the same as the MB-pol intramolecular component. For comparison, we also computed the RDFs using the MB-pol and q-tip4p/f force fields at the same condition.