Bayesian Modeling of Polarizable Water: Lessons for Force Field Development

Computer simulation of molecular dynamics is useful only when the accuracy of models and the uncertainty of model predictions can be quantified. The accuracy of simulations can usually be established by comparing with reference data from experiments or quantum chemistry. However, trade-offs in designing simulation models often lead to property-dependent accuracy. For instance, models trained to reproduce liquid-phase density and enthalpy of vaporization do not automatically reproduce solid or gas-phase properties. This means that the method chosen for force field design introduces uncertainties in the final model, which are often ignored. Although the statistical uncertainty due to stochastic simulations is usually accounted for, the uncertainty due to parameters is harder to quantify. In a recent paper( npj Comput. Mater. 11 (2025), 366)., we developed Bayesian three-point water models based on synthetic likelihoods. The results were shown to depend strongly on the reference observables used. Here, we extend this work to polarizable models by applying Bayesian inference to the SWM4-NDP model ( Chem. Phys. Lett. 418 (2006), 245). We consider two van der Waals functional forms, and subject their corresponding parameters to Markov-chain Monte Carlo sampling. The two functional forms have similar performance despite the differences in their complexity, but favor different inference observables. In addition to this, we find that explicit polarization does not significantly reduce parameter uncertainty for the chosen observable set. We then compare models identified by Bayesian inference with models in which the van der Waals potential was trained on gas-phase dimer energies from symmetry-adapted perturbation theory using the Alexandria Chemistry Toolkit ( Digit. Discovery 4 (2025), 1925).Intriguingly, the Lennard-Jones 12-6 parameters trained on the gas-phase dimers are very similar to the optimal parameter set from the Bayesian inference, and the resulting models have similar accuracy. Replacing the Lennard-Jones 12-6 potential with the Wang-Buckingham potential ( J. Chem. Theory Comput. 9 (2013), 452)and training on gas-phase data yields an unstable model, whereas the corresponding optimal Bayesian model gives reasonable results. A further analysis shows that observables respond in different ways to changes in force field parameters, highlighting the need to consider a large number of observables during force field training. Implications for model development as well as uncertainty quantification are discussed.