The water molecule’s electronic Cartesian multipole moment and polarizability tensors have been fitted with Gaussian process regression to the internal coordinates and are used to evaluate accurate electrostatic, induction, and dispersion energy components between flexible molecules. The model yields a handful of damping and scaling parameters that were adjusted for the energy components to agree with 2-body symmetry-adapted perturbation theory decomposition and then fine-tuned in order for the total energy to agree with CCSD(T) for small clusters. We present a simple algorithm for rotating symmetric Cartesian tensors and employ a dispersion potential based on multipole polarizabilities. At short range, the 2- and 3-body potential energy was corrected to CCSD(T) accuracy using Gaussian approximation potentials. The radial distribution function and self-diffusion coefficient obtained with molecular dynamics simulations agree well with experiments.