here we introduce the analytical 
 equations that yield the positions 
 and values of the three point charges that best reproduce the three 
 lowest order multipole moments of the water molecule .
the lowest three 
 nonzero multipole moments of the water molecule are the dipole that 
 is represented by one independent component (  ) , the quadrupole 
 defined by two independent components ( q0 , qt ) , and the octupole defined by two 
 independent components ( 0 , t ) . in the coordinate system shown in figure 2 , these moments are related to the cartesian components 
 of the traceless multipole moments of water molecule as  = 
 z , q0 = qzz , qt = 1/2(qyy  qxx ) , 0 = ozzz , and 
 t = 1/2(oyyz  oxxz ) ( see si ) .
the optimal point charges are calculated so 
 that these moments 
 are sequentially reproduced , starting with the lowest order moments .
the dipole and the quadrupole moments are reproduced 
 exactly by requiring123where z2 , z1 , y and q are the independent 
 unknown parameters that characterize the three 
 point charge model ( see figure 2 ) .
the above 
 set of equations is solved to find three geometrical parameters of 
 the water model ( z2 , z1 , and y).45 this leaves only 
 one unknown parameter , the charge value q , which 
 we calculate by using two additional equations that relate the charge 
 distribution parameters to the octupole moment components so that 
 the octupole moment is optimally reproduced ( see si ) .
the calculations of thermodynamic 
 and dynamical bulk properties 
 were done based on standard equations in the literature ( see si for details ) . unless specified otherwise , 
 we use the following md simulations protocol . simulations in the npt 
 ensemble ( 1 bar , 298.16 k )
a cubic box with edge length of 30 
  was filled with 804 water molecules .
long - range electrostatic interactions , calculated via the 
 particle mesh ewald ( pme ) summation , and the van der waals interactions 
 were cut off at distance 8  .
md simulations were conducted with 
 a 2 fs time step ; all intramolecular geometries were constrained with 
 shake .
the npt simulations were performed using langevin thermostat 
 with coupling constant  = 2.0 ps , and a 
 berendsen barostat with coupling constant of 1.0 ps for equilibration and 3.0 ps for production . 

the duration of production runs vary between 1 to 65 ns , 
 depending on the properties ( see si ) . to mitigate uncertainties due to conformational variability , the 
 20 test molecule were randomly selected from a subset of 248 highly 
 rigid molecules .
explicit solvent free 
 energies calculations ( via thermodynamic integration ) were performed 
 in gromacs 4.6.5 using the gaff small molecule parameters ( see si for further details ) .
the predictive power of models 
 against experimental data was validated 
 using a scoring system developed by vega et al . for a calculated property x and a corresponding 
 experimental value of xexp ,
the assigned 
 score is obtained as6where 
 the tolerance ( tol ) is assigned to 0.5% 
 for density , position of the first peak of the rdf and for heat of 
 vaporization , 5% for height of the first peak of the rdf , and 2.5% 
 for the remaining properties .
the quality score assigned to each test 
 model is equal to the average of the scores in bulk properties considered .
 simplified 
 classical water models are currently an indispensable 
 component in practical atomistic simulations . yet , despite several 
 decades of intense research , these models are still far from perfect . 
 
 presented here is an alternative approach to constructing widely used 
 point charge water models . 
 in contrast to the conventional approach , 
 we do not impose any geometry constraints on the model other than 
 the symmetry . 
 instead , we optimize the distribution of point charges 
 to best describe the  electrostatics  of the water molecule . 
 the resulting  optimal  
 3-charge , 4-point rigid water 
 model ( opc ) reproduces a comprehensive set of bulk properties significantly 
 more accurately than commonly used rigid models : average error relative 
 to experiment is 0.76% . 
 close agreement with experiment holds over 
 a wide range of temperatures . 
 the improvements in the proposed model 
 extend beyond bulk properties : compared to common rigid models , predicted 
 hydration free energies of small molecules using opc are uniformly 
 closer to experiment , with root - mean - square error < 1 kcal / mol . 
Computational Methods